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Foreword 


This  research  program  has  focused  on  enhancing  the  understanding  of  the  internal  flows  in  injector  passages 
and  their  contributions  to  the  subsequent  spray  formation  outside  the  orifice.  During  the  past  three  years, 
progress  has  been  made  in  the  following  aspects.  First,  the  investigation  of  the  unsteadiness  in  orifice 
massflow  caused  by  the  hydrodynamic  instability  of  the  vena-contract  a  or  the  presence  of  cavitation  in  this 
region  through  a  parametric  study  of  an  axisymmetric  orifice.  Second,  the  k  —  to  turbulence  model  and  the 
homogeneous  pseudo  density  model  has  been  combined  to  approximate  two-dimensional  and  axisymmetric 
turbulent  cavitating  flows.  Although,  the  turbulence  model  generate  a  steady  attached  cavity  at  the  inlet 
corner,  it  shows  an  improvement  in  the  prediction  of  mass  flow  through  an  orifice.  Finally,  three  dimensional 
laminar  calculations  were  conducted  to  provide  a  better  understanding  of  cross-flow  effects.  The  effect  of 
cross  flow  velocity  on  orifice  mass  flow  is  investigated  and  the  three  dimensional  model  predicts  a  satisfactory 
results  compared  with  experimental  measurements. 
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Figure  1:  Computational  domain  and  mesh  for  orifice  flow. 


1  Research  Objectives 

Notable  changes  in  spray  structure  are  evident  under  the  variations  of  the  cavitation  extent  inside  the  injector 
passage  [^2,3]  fj^nce,  the  understanding  of  the  complex  cavitation  phenomena  present  in  an  injector  passage 
is  fundamental  to  the  subsequent  jet  atomization.  For  this  purpose,  a  focused  research  effort  has  been 
conducted  to  develop  models  capable  of  providing  quantitative  information  regarding  the  cavitation  process. 
The  models  have  centered  on  the  use  of  Marker  and  Cell  finite  volume  method  as  a  means  to  provide  accurate 
description  of  the  complex,  and  arbitrary  unsteady  conditions. 

The  homogeneous  pseudo  density  model,  developed  by  Chen  and  Heister  which  describes  the  non¬ 
equilibrium  interaction  between  liquid  and  bubbles  provides  the  basis  of  this  research.  It  assumes  two  phases 
to  be  fully  mixed  on  the  sub-grid  level  and  the  local  density  is  a  measure  of  void  fraction.  This  simplification 
circumvents  the  great  challenge  of  tracing  each  individual  bubble  which  is  implausible  at  current  stage,  but 
instead  focuses  on  the  global/integral  dynamics  of  bubble  clusters.  Based  on  this  pseudo  density  formulation, 
a  series  of  CFD  program  have  been  developed  with  increasing  physics  and  complexity:  (1)  two-dimensional 
and  axisymmetric  laminar  (2)  two-dimensional  and  axisymmetric  turbulent^^^ ,  (3)  three-dimensional 
laminar^^l  These  CFD  programs  were  used  to  predict  the  evolution  of  cavitation  and  its  effect  on  orifice 
mass  flow. 


2  Summary  of  The  Most  Important  Results 

Three  major  tasks  have  been  accomplished  during  the  three  years  research  program.  The  unsteady  effect  of 
the  instability  of  the  vena-contracta  or  cavitation  in  this  region  on  massflow  is  addressed  in  Section  3.1.  The 
influence  of  cavitation  on  flow  turbulence,  and  turbulent  orifice/slot  internal  flow  is  described  in  Section  3.2. 
Finally,  the  results  from  an  unsteady  three-dimensional  two-phase  model  are  provided  in  Section  3.3. 

Figure  1  shows  a  typical  mesh  for  an  axisymmetric  orifice  or  a  two-dimensional  slot.  Grid  refinement 
study  were  performed  for  both  laminar  and  turbulent  calculations.  It  is  verified  that  a  mesh  employing  140 
grid  points  in  the  axial  direction  and  60  points  in  the  radial  direction  is  adequate  to  resolve  salient  flow-field 
structures  for  laminar  calculation  .  For  turbulent  calculation  the  number  of  grid  points  in  radial  direction 
is  increased  to  90.  Constant  pressure  boundary  conditions  are  imposed  on  inflow  and  outflow  boundaries; 
no-slip  conditions  are  imposed  along  walls,  and  symmetry  conditions  are  imposed  along  the  centerline.  To 
approximate  the  velocity  at  the  inflow  boundary,  we  employ  a  sink  at  the  origin.  The  strength  of  the  sink  is 
updated  during  each  time  step  by  the  conservation  of  mass  flow  rate  through  the  nozzle  passage. 

2.1  Laminar  Parametric  Study  on  an  Axisymmetric  Orifice 

A  series  laminar  simulations  were  conducted  to  assess  the  unsteady  flow  perturbations  brought  about  by 
the  vena-contracta  in  single  and  two-phase  regimes.  The  influence  of  Reynolds  number,  supply  pressure, 
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discharge  pressure,  orifice  L/D  and  inlet  rounding  have  been  investigated.  In  this  report,  the  unsteadiness 
in  orifice  massflow  and  the  effect  of  inlet  rounding  are  described. 

2.1.1  Unsteadiness  in  Massflow  and  Cavitation  Length 

The  predicted  orifice  massflow  characteristics  shows  a  quasi-periodic  oscillation  under  both  cavitating  and 
non-cavitating  conditions.  Under  the  presence  of  cavitation,  the  oscillation  of  the  extent  of  cavitation  region 
is  the  primary  reason  that  causes  the  unsteadiness  of  the  massflow  through  orifice.  Figure  2  shows  the  time 
history  of  cavitation  length  {Lc)  and  discharge  coefficient  {Cd)  for  a  typical  cavitating  conditions. 


Figure  2;  Typical  quasi-periodic  behavior  of  cavitation  length  Lc  and  orifice  discharge  coefficient  Cd 

By  performing  a  Fourier  transform  on  the  Lc  and  Cd  histories,  one  can  obtain  the  principle  frequency  of 
oscillation  of  each  of  these  parameters,  fc  and  /d  •  We  should  point  out  that  the  Fourier  transforms  typically 
show  energy  in  2-3  discrete  frequencies^^’^®^.  The  bulk  of  the  energy  is  contained  in  the  primary  harmonic, 
other  frequencies  which  appear  to  be  higher  harmonics  tend  to  contain  much  less  energy.  In  addition,  the 
oscillation  of  Lc  and  Cd  is  coupled  tightly.  This  point  is  evidenced  in  figure  3  which  shows  the  change  in  fc 
and  /d  with  variations  in  K  at  fixed  Reynolds  number  for  an  orifice  with  L/D  S.  The  point  of  cavitation 
inception  is  shown  on  the  figure;  the  region  to  the  left  of  this  point  is  where  cavitation  is  evident  in  the  flow 
field.  Note  in  the  plot,  fc  and  /d  are  almost  identical  under  the  presence  of  cavitation. 

It  may  be  possible  to  make  use  of  this  result  to  design  passive  oscillations  of  a  desired  frequency.  In 
general,  the  frequencies  are  consistent  with  the  time  it  takes  a  fluid  element  to  traverse  the  length  of  the 
orifice  passage.  Atomizers  which  produce  fine  sprays  tend  to  produce  droplets  at  much  higher  frequencies 
than  those  noted  in  figure  3.  Partial  cavitation  may  be  used  to  “pump”  instabilities  for  atomizers  operating 
at  more  modest  pressure  drops  such  as  a  flow  produced  in  inkjet  printing  applications. 

2.1.2  Effect  of  Inlet  Rounding 

Several  experiments  and  numerical  simulationsf^^  have  noted  that  inlet  rounding  can  delay  the  occurrence  of 
cavitation  and  increase  discharge  coefficient  through  reductions  in  contraction  losses.  However,  the  effects 
of  rounding  on  the  unsteady  behavior  has  not  been  investigated  in  any  great  detail.  For  this  reason,  the 
unsteady  characteristics  have  been  studied  assuming  the  inlet  is  rounded  with  a  given  radius  circle.  Figure  4 
depicts  the  level  of  massflow  variations  at  various  cavitation  numbers  for  three  levels  of  inlet  rounding 
{r/R  =  0,  0.05,  0.1)  and  the  change  of  average  discharge  coefficient  vs.  inlet  rounding  at  fixed  cavitation 
number  K  =  1.56.  The  cavitation  onset  point  is  noted  on  the  left  plot;  cavitating  results  lie  to  the  left  of 
this  point  since  decreasing  K  corresponds  to  increased  supply  pressure.  Results  in  the  figures  show  that 
massflow  variations  are  reduced  dramatically  for  even  a  small  amount  of  inlet  rounding,  while  the  average 
mass  flow  rate  is  increased. 
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Figure  3:  Frequency  of  discharge  coefficient  (/p)  and  cavitation  length  {fc)  vs.  cavitation  number;  L/D  = 
8,  P2  =  1  atm. 


Figure  4:  The  effect  of  inlet  rounding  on  the  fluctuation  of  discharge  coefficient  ACd  and  its  mean  value  C d 

2.2  Turbulent  Injector  Internal  Flows 

A  few  axisymmetric  turbulent  simulations  were  made  and  the  predicted  discharge  coefficient  is  compared 
with  the  measurements  by  Nurick^^^^.  The  two  dimensional  code  was  used  to  assess  the  turbulent  cavitating 
flow  through  a  slot.  The  computed  cavitation  extension  and  turbulent  quantities  is  compared  with  available 
experimental  measurements.  The  results  are  described  in  the  following  sections. 

2.2.1  Discharge  Coefficient  Prediction 

Figure  5  shows  the  results  of  the  discharge  coefficient  comparison  on  an  orifice  of  L/D  =:  6,Z)  =  3.18mm 
under  a  back  pressure  of  P2  =  13.8ps2.  In  the  figure,  “Laml”  represents  a  laminar  calculation  without 
sink-inflow  velocity  boundary  condition,  “Lam2”  represents  a  laminar  calculation  with  an  approximated 
inflow  velocity  from  an  artificial  sink  as  described  in  the  previous  section,  and  the  line  denoted  as  “Turb” 
is  calculated  from  the  turbulence  model.  As  shown  in  the  plot,  the  sink-inflow  velocity  treatment  greatly 
improves  the  prediction  of  discharge  coefficient.  There  is  an  increase  in  the  magnitude  oi  Cd  of  4.7%  over 
the  zero  inflow  velocity  treatment.  The  turbulence  model  shows  a  further  improvement  and  gives  results 
somewhat  closer  to  the  experimental  data.  The  differences  between  the  turbulence  model  and  laminar 
predictions  on  Cd  might  be  explained  by  the  exit  velocity  profiles  from  both  types  calculation.  The  turbulent 
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Figure  5:  Discharge  coefficient  Cd  comparison  with  experimental  results;  L/D  =  6,  D  =  3.18mm;  and  the 
comparison  of  velocity  profile  at  the  exit  for  laminar  and  turbulent  solutions 


velocity  profile  is  fuller  than  the  laminar  one.  Due  to  this  feature,  the  turbulence  model  yields  a  larger 
Cd  by  1.6%.  In  the  calculations  a  sharp-edged  inlet  is  assumed,  however  the  amount  of  inlet  rounding 
on  the  experimental  hardware  is  not  known.  As  indicated  in  figure  4  small  degree  of  inlet  rounding  can 
lead  to  substantial  increase  in  discharge  coefficient.  It  is  estimated  that  an  inlet  rounding  less  than  0.025 
would  be  sufficient  for  the  laminar  model  to  generate  a  discharge  coefficient  that  is  high  enough  to  match 
the  measurement.  This  factor  might  be  the  primary  reason  for  the  model’s  consistent  underprediction  of 
discharge  coefficient. 


Figure  6:  Overlays  of  numerical  results  and  photographic  time  shot  of  cavitation  region,  (AP  36psi,  LfD  = 
10.714).  Here,  p  —  0.98  on  the  outermost  contour  with  gradients  of  0.2  between  contours 

2.2.2  Cavitation  Extent  Comparison 

The  2-D  code  was  used  to  simulate  the  flow  through  slots  with  the  geometry  used  by  Henry^^^  and  Sanchez^^^^. 
Although  the  turbulent  code  produces  a  constant  value  of  the  cavitation  length,  the  turbulent  calculation 
generates  a  single  contiguous  cavitation  region  near  the  inlet  corner,  consistent  with  experimental  observation. 
Figure  6  shows  density  contours,  which  denote  the  cavity  region,  obtained  from  laminar  (the  upper  one) 
and  turbulent  (the  lower  one)  calculations,  overlaying  one  photographic  snapshot  by  Henry Although  the 
laminar  simulation  results  in  an  overall  cavitation  extent  consistent  with  experiment,  it  indicates  two  separate 
regions  of  cavitation.  The  turbulence  model  improves  on  this  point  by  generating  a  single  cavitation  region 
which  appears  to  be  quite  consistent  with  experimental  results  both  in  axial  and  cross-stream  directions. 
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2.2.3  Turbulent  Kinetic  Energy  (TKE) 

The  turbulence  model  predicts  highest  TKE  at  the  rear  end  of  cavity  region,  the  wake  of  the  cavity,  for  all 
operating  conditions.  It  is  known  that  the  wake  where  bubble  detachment  and  collapses  accompanied  with 
pressure  recovery  occurs  is  the  most  dynamic  region.  The  computed  turbulence  quantity  is  compared  with 
the  measurement  of  Ruiz  and  Figure  7  shows  turbulence  intensities  in  the  streamwise  and  normal 

direction  and  Reynolds  shear  stress  as  well  as  the  experimental  data  at  different  streamwise  locations. 

Here,  R^x  =  pu  u”  and  Ryy  =  ^ pv*  v” ,  representing  turbulence  intensity  in  x  and  y  directions  are  the 
dimensionless  mass  weighted  velocity  fluctuation  in  streamwise  and  normal  directions  ,  and  Rxy  —  is 

the  mass  averaged  Reynolds  shear  stress.  The  cavitation  extent  is  also  shown  in  this  plot. 

1  1 

Near  the  rear  end  of  cavity  region  the  calculation  predicts  Rix  yRyy  and  Rxy  are  slightly  larger  under 
cavitating  conditions  than  under  non-cavitating  conditions  due  to  flow  reattachment.  It  should  be  mentioned 
the  k  —  oj  model  does  not  include  any  production  terms  due  to  bubble  collapse.  The  only  mechanism  to 
capture  turbulence  generation  for  the  turbulence  model  is  the  occurrence  of  a  strong  strain  rate.  However, 
due  to  the  presence  of  cavitation  the  pseudo  density  model  produces  a  larger  strain  rate  flow  field.  It  has 
been  observed  that  the  collapse  of  bubbles  is  a  source  of  vortex  generation^^^^ . 


2.3  The  Effect  of  Cross  Flow 

The  fully  three-dimensional  two-phase  Niavier-Stokes  solver  is  utilized  to  simulate  an  orifice  flow  driven  by 
the  cross  flow  in  a  manifold.  Figure  8  shows  the  schematic  representation  of  the  three-dimensional  manifold 
cross  flow.  In  the  figure,  Vi  is  the  cross  flow  velocity,  and  V2  is  ideal  discharge  velocity  at  the  exit,  calculated 
from  Bernoulli’s  equation.  The  effect  of  the  cross  flow  velocity  on  cavitation  length  and  discharge  coefficient 
are  investigated.  The  computed  discharge  coefficient  is  compared  with  Strakey  and  Talley measurement. 

2.3.1  Discharge  Coefficient 

Figure  9  shows  the  comparison  of  the  computed  average  discharge  coefficient  and  the  experimentally  measured 
discharge  coefficient.  Note  the  average  value  is  obtained  by  sampling  the  data  after  it  reaches  a  steady 
oscillation  state.  Although  the  number  of  computed  data  is  fewer  than  that  of  the  experimental  results 
because  of  the  high  expense  of  the  computation,  two  patterns  can  still  be  observed  from  the  plot.  First,  for 
a  given  cross  flow  velocity,  as  cavitation  number  increases  the  discharge  coefficient  increases  at  first.  After 
reaching  the  maximum  value,  discharge  coefficient  then  begins  to  drop  off.  The  maximum  value  in  discharge 
coefficient  occurs  when  the  cavitation  number  is  about  1.8,  a  value  close  to  the  experimentally  measured 
inception  index.  The  reason  for  that  is  the  following.  The  mass  flow  rate  is  proportional  to  the  pressure 
drop  through  the  orifice  (increasing  with  decreasing  cavitation  number)  before  the  inception  of  cavitation. 
When  the  pressure  loss  is  large  enough  to  cause  cavitation  to  occur,  the  mass  flow  rate  begins  to  drop  due  to 
the  slipstream  effect  of  the  cavity  as  pointed  out  by  Bunnell  and  Heister^^®^ .  The  more  significant  cavitation 
becomes,  the  lower  the  mass  flow  rate  is.  The  limited  computed  data  samples  reproduce  this  feature  as  the 
experimental  counterparts  do,  except  that  they  overpredict  the  discharge  coefficient  in  general,  especially 
when  Vi  =  S,9m/s  and  K  ~  6.0.  The  overprediction  might  be  due  to  the  lack  of  turbulence  in  the  modeling 
or  because  the  grid  is  not  fine  enough^^^ 

The  second  trend  observed  in  figure  9  is  with  regard  to  the  effect  of  cross  flow  velocity.  Both  the 
computation  and  experiment  show  that  as  cross  velocity  increases,  discharge  coefficient  decreases.  This 
phenomenon  can  be  attributed  to  the  increase  in  velocity  ratio,  the  ratio  of  the  cross  flow  velocity  to  the 
ideal  orifice  discharge  velocity  With  increasing  velocity  ratio,  less  amount  of  fluid  passes  through 

the  orifice. 

2.3.2  Internal  Flow  Structure 

Soteriou  et  al.^^^1  have  found  strong  vortex  structures  exist  inside  the  orifice  when  cavitation  occurs  based 
on  their  observation  on  a  large  scale  nozzle.  Vortex  interaction  can  also  be  observed  under  non-cavitating 
conditions.  The  vortices  inside  the  nozzle  sometimes  intertwine  with  each  other.  Figure  10  shows  the 
photographic  shot  of  a  vortex  intertwining  structure  inside  an  orifice,  and  the  vortex  structure  of  a  diesel 
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Figure  7:  Turbulence  intensities  and  mass- averaged  Reynolds  stress  comparison  with  experimental  results 
at  different  axial  locations,  with  and  without  cavitation  (turbulent  calculation) 
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Figure  8:  3-D  manifold  cross  flow  model  and  key  parameters 
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Figure  9:  The  effect  of  cross  flow  velocity  on  discharge  coefficient  and  the  comparison  with  experiment 
measurement 


spray.  Contrary  to  the  traditional  theory  that  the  aerodynamic  interfacial  shear  is  the  dominant  force  that 
causes  the  breakup  of  the  fuel  injected  from  the  orifice,  Soteriou  et  al.  indicated  that  the  primary  factors 
inducing  atomization  are  associated  with  the  characteristics  of  the  internal  flow  inside  the  nozzle. 

Figure  11  shows  the  various  computed  vortex  lines  inside  the  orifice  at  one  instant  in  time  for  the  condi¬ 
tions  with  and  without  cavitation.  At  other  instants,  both  flow  conditions  show  a  vortex  structure  similar 
to  what  are  shown  in  these  two  figures.  As  seen  from  the  figure,  the  difference  between  the  vortex  structure 
under  cavitating  condition  and  that  under  non-cavitating  condition  is  significant.  For  non-cavitating  case, 
in  the  front  section  of  the  orifice,  the  cross  or  vertical  components  of  the  vorticity  are  dominant.  The  few 
vortex  lines  shown  in  the  figure  are  nearly  parallel  to  each  other  and  perpendicular  to  streamwise  direction. 
Based  on  these  representing  vortex  lines,  we  believe  vortex  sheet-like  structures  exist  near  the  wall  region  in 
the  front  section  of  the  orifice.  In  the  rear  section  of  orifice,  the  streamwise  component  of  the  vorticity  start 
to  be  dominant.  The  vortex  sheet  starts  to  roll  up  and  several  vertex  lines  intertwine  and  form  a  strong 
vortex  structure  on  the  windward  side  of  orifice.  For  the  flow  with  cavitation,  vortex  interaction  begins 
earlier  on  the  leeward  side  and  two  primary  vortex  structures  are  formed  at  the  exit  of  orifice.  These  vortex 
structures  indicates  that  swirls  with  rotating  axis  in  line  with  streamwise  direction  exist  at  the  exit. 

The  role  played  by  the  internal  vortex  structure  on  spray  atomization  is  not  well  understood  since  it 
is  both  a  numerical  and  an  experimental  new  discovery.  Recently  some  authors^^’^’^^^  have  suggested  that 
cavitation  inside  the  orifice  is  an  important  factor  that  promotes  atomization.  This  might  be  attributed  to 
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Figure  10:  Experimental  snapshots  of  vortex  structure  inside  an  orifice  (on  the  left)  and  of  diesel  sprays  (on 
the  right).  [From  Soteriou  et  used  by  permission] 


Figure  11:  Vortex  structure  inside  the  orifice  under  cavitating  {K  =  1.2  on  the  left)  and  non-cavitating 
(/\  =  6.0,  on  the  right)  conditions,  Vi  =  6.0m/5.  Green  lines  are  vortex  lines  outlined  parallel  to  local 
vorticity  vector. 


the  difference  made  by  cavitation  on  the  internal  flow  structure  and  the  magnitude  of  the  vorticity  vector 
at  the  exit  plane.  The  simulation  shows  that  there  is  a  stronger  swirl  on  the  leeward  side  for  cavitating 
flow.  Once  the  fuel  flows  out  of  the  orifice,  without  the  constraint  of  the  wall  boundary  the  stronger  velocity 
components  in  cross  and  vertical  directions  tend  to  cause  the  jet  to  break  up  earlier. 


3  Technology  Transfer 

We  continue  to  work  closely  with  Cummins  Engines  in  transfering  results  from  this  research  to  the  field. 
In  late  1999,  we  completed  a  series  of  fully  3-D  unsteady  calculations  of  one  of  Cummins  injector  designs. 
These  demanding  calculations  were  not  totally  coupled  to  the  inflow  from  the  injector  plunger  and  we  have 
plans  to  continue  the  work  to  include  this  coupling.  We  have  also  had  several  conversations  with  officials  at 
Detroit  Diesel  and  Catapiller  at  national  meetings  and  continue  to  inform  the  community  of  our  findings. 
Politically,  it  is  difficult  for  us  to  work  directly  with  these  other  organizations  due  to  conflict  of  interest 
concerns  by  Cummins. 
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Abstract 

A  series  of  parametric  simulations  have  been  conducted  to  investigate  the  unsteady  internal  flow  behavior 
in  plain  orifice,  or  “pressure”  atomizers.  The  unsteady  behavior  is  attributed  to  instabilities  in  the  vena 
contracta  and  the  presence  of  cavitation  in  this  region.  Even  though  the  present  calculations  are  laminar, 
they  exhibit  turbulent-like  characteristics  due  to  this  instability.  Orifice  massflow  is  shown  to  be  periodic  in 
most  instances  with  fluctuations  occurring  at  frequencies  consistent  with  the  orifice  transit  time.  Cavitation 
is  shown  to  enhance  fluctuations;  in  some  cases  quite  dramatically.  Results  are  presented  for  a  range  of  flow 
conditions  and  orifice  geometries. 

Introduction 

The  plain-orifice,  or  pressure  atomizer  represents  the  simplest  solution  for  atomizing  liquids  in  many  ap¬ 
plications.  By  simply  drilling  a  small  hole  in  a  piece  of  material,  adequate  atomization  can  be  obtained  in 
many  instances  for  fluids  of  low  or  modest  viscosity.  For  this  reason,  the  pressure  atomizer  has  very  wide 
usage  and  has  been  the  subject  of  research  studies  for  well  over  100  years ^ .  Innumerable  works  have  focused 
on  the  droplet /spray  produced  by  these  devices  and  modern  instrumentation  provides  accurate  measures  of 
droplet  sizes  and  velocities  in  dilute  regions  downstream  from  the  orifice  exit. 

It  has  long  been  recognized  the  design  of  the  orifice  flow  passage  (both  length  and  diameter)  has  important 
implications  on  the  type  of  spray  pattern  produced  by  the  orifice^’^.  Through  changes  in  orifice  design,  a  high 
pressure  drop  atomizer  can  be  manufactured  to  produce  a  very  fine  spray  or  a  water-jet  cutter.  The  tendency 
for  the  flow  to  cavitate  on  the  inner  orifice  lip  has  also  been  inferred  as  a  basic  atomization  mechanism  by 
some  authors^“^.  These  factors  have  led  numerous  researchers  to  begin  to  study  both  experimentally^" 
and  analytically^ the  flow  pattern  inside  the  orifice  passage  in  hopes  to  gain  insight  on  its  effect  on  spray 
attributes. 

Motivated  primarily  by  the  diesel  injector  application,  recent  efforts  have  focused  on  the  impact  of  cavita¬ 
tion  on  internal  flow  and  spray  characteristics.  Figure  1  highlights  the  simple  geometry  for  a  typical  pressure 
atomizer  with  upstream  and  downstream  pressures,  Pi  and  P2,  respectively.  The  large  fluid  acceleration 
near  the  inlet  lip  leads  to  a  locally  low  pressure  region;  a  vena-contracta  is  formed  and  cavitation  will  first 
appear  in  this  area.  Numerous  experiments  have  shown  the  cavitation  region  to  be  inherently  unsteady;  the 
collapsing  bubbles  at  the  aft  end  to  the  region  cause  local  pressure  rises  which  are  fed  upstream  and  effect 
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subsequent  shapes  of  the  cavity.  Recent  experiments^^  have  shown  that  this  “partial  cavitation”  behavior 
enhances  instabilities  in  the  jet  produced  by  the  orifice  thereby  promoting  atomization. 

Even  under  non-cavitated  conditions,  the  vena-contracta  is  still  present  and  is  subject  to  instabilities 
arising  from  the  abrupt  pressure  rise  at  the  end  of  the  recirculation  zone.  As  the  pressure  drop  Pi  “  P2 
increases,  the  flow  in  the  passage  will  eventually  reach  the  hydraulic  flip  condition  in  which  the  vena  contracta 
extends  the  entire  length  of  the  orifice.  Notable  changes  in  spray  structure  are  evident  under  these  conditions; 
atomization  is  generally  much  poorer  than  under  partial  cavitation  conditions. 

Numerical  modeling  of  these  orifice  flows  is  challenging  because  one  must  normally  account  for  unsteadi¬ 
ness  and  the  possibility  of  cavitation  within  the  passage.  Recently,  several  authors^^d^, 19,22  approached 

this  problem  with  a  homogeneous,  or  “pseudo-density”  formulation  in  which  the  single  phase  Navier  Stokes 
equations  are  solved  on  a  fixed  computational  mesh.  An  additional  constitutive  relation  is  required  for  the 
pseudo-density  in  these  schemes.  These  schemes  are  powerful  in  that  they  can  normally  provide  economical 
calculations  of  flows  in  which  there  are  simply  too  many  bubbles  to  track  individually. 

In  this  paper,  we  apply  a  recently  developed  homogeneous  flow  model  to  assess  the  unsteady  flow  per¬ 
turbations  brought  about  by  the  vena-contracta  in  single  and  two-phase  flow  regimes.  The  following  section 
provides  a  brief  description  of  the  model  and  its  validation  against  experimental  data.  A  large  number  of 
parametric  simulations  are  then  summarized,  followed  by  conclusion  from  the  studies. 


Modeling  Description 


Since  the  numerical  model  has  been  described  in  great  detail  in  other  works^^”^^,  we  will  provide  just 
a  brief  overview  here.  We  assume  a  laminar,  axisymmetric,  incompressible  flow  and  solve  the  Navier- 
Stokes  equations  on  a  fixed  structured  mesh  using  the  Marker  and  Cell  algorithm.  Liquid  density  (p/), 
the  “Bernoulli”  velocity  in  the  orifice  {v  =  >/2(Pi  ~  P2)/pi),  and  the  orifice  diameter  (D)  are  chosen  as 
dimensions.  In  this  case,  the  two  dimensionless  parameters  are  the  Reynolds  number  {Re)  and  cavitation 
number  (A"): 


Re  — 


pivD 


K  = 


Pl-Pv 


(1) 


pi  Pi  “  P2 

where  Py  is  the  fluid  vapor  pressure  and  pi  is  the  liquid  viscosity.  The  viscosity  in  the  two-phase  mixture  is 
generally  based  on  the  local  void  fraction^®,  a: 


p  =  apg  -i-  (1  -  a)pi 


(2) 


Since  pg  «  Pi  and  Pg  «  pi  for  most  liquid/vapor  mixtures,  we  can  effectively  neglect  the  influence  of 
gas  phase  viscosity.  Under  these  assumptions,  the  pseudo-density  of  the  mixture  is  p  ==  1  —  a  and  Eq.  2 
becomes: 

n  =  PHI  (3) 

Figure  2  shows  a  typical  mesh  for  an  orifice  L/D=6.  The  mesh  employed  140  grid  points  in  the  axial 
direction  and  60  points  in  the  radial  direction.  Recent  grid  convergence  studies  have  verified  that  this  mesh 
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is  adequate  to  resolve  salient  flowfield  structures'^  for  non-cavitating  or  cavitating  flows.  Constant  pressure 
boundary  conditions  are  imposed  on  inflow  and  outflow  boundaries;  no-slip  conditions  are  used  along  walls, 
and  symmetry  conditions  are  imposed  along  the  orifice  centerline. 

The  two-phase  treatment  requires  a  constitutive  relation  for  the  pseudo-density,  p,  in  order  to  obtain  clo¬ 
sure  for  the  governing  equations.  The  current  treatment^^  is  based  on  the  dynamic  response  of  a  monodisperse 
bubble  field  to  changes  in  pressure  and  inertial  forces  per  the  Rayleigh- Plesset  equation.  The  treatment  ne¬ 
glects  surface  tension,  bubble  coalescence,  bubble  breakup,  and  slip  between  phases.  However,  the  approach 
does  apply  for  arbitrary  void  fraction  (or  pseudo-density)  and  we  do  not  assume  a  disperse  bubble  field  used 
by  many  researchers.  These  assumptions  greatly  simplify  the  physics  of  the  flow,  but  still  provide  a  gross 
response  to  the  spatial  and  time  varying  pressure  field.  The  resulting  constitutive  relation  takes  on  a  form 
similar  to  the  Rayleigh- Plesset  equation  in  that  changes  in  pseudo-density  are  governed  by  pressure  and 
inertial  forces: 

-py  ^  +  a  -  P  )+ 

-  a'2  -  a' - 1  l+Aa'+a^  ('4'j 

'•  6a'3(l-a'3)  6a'2(2  +  a')(H-a' 

Here,  a'  —  -^1  —  p,  and  P  is  the  local  pressure  returned  by  the  Navier-Stokes  solver.  Finally,  in  Eq.  4 
Lq  is  a  non-dimensional  characteristic  length  scale: 


Here,  Uo  is  the  non-dimensional  site  density.  Letting  the  dimensional  site  density  be  represented  by  Uoj  then 

Uo  =  rioD^  (6) 

The  value  of  Uo  is  a  subject  of  debate  since  it  is  currently  impossible  to  measure  submicron  scale  bubbles 
which  can  presumably  serve  as  nucleations  sites.  Current  estimates  lie  in  the  range  of  10®  -  10^^ sites/ np 
for  small  scale  internal  flows.  Recent  parametric  analyses^ ^  indicate  that  results  are  not  sensitive  to  this 
parameter;  a  value  of  10^^  sites/m®  was  used  for  all  calculations  presented  herein. 

As  the  number  of  sites  decrease  with  the  volume  of  liquid,  Lq  grows  and  the  pressure  difference  term  in 
Eq.  4  becomes  less  dominant.  This  feature  accounts  for  hydrodynamic  non-equilibrium  effects  in  which  the 
inertia  of  the  liquid  (second  term  on  RHS  of  Eq.  4)  becomes  significant.  It  is  this  feature  that  accounts  for 
scaling  effects  which  are  well  established  in  cavitating  flows. 

Model  Validation 

The  model  has  had  extensive  validation  on  a  variety  of  cavitating  flows.  Cavitation  extent  {Lc  in  Fig.  1) 
matches  experimental  results  for  flows  over  axisymmetric  headforms^®^^^  In  addition,  extensive  comparisons 
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have  been  made  on  experiments  conducted  in  a  high  aspect  ratio  cavitating  slot  flow  using  a  2-D  approx¬ 
imation  for  the  flowfield.  In  this  case,  cavitation  lengths^d2  compare  well  with  measurements  from  Dr. 
Collicott’s  research  group  at  Purdue. 

For  the  present  studies,  comparisons  were  made  with  cavitation  inception  pressure  measurements  of 
Bergwerk^,  with  discharge  coefficient  measurements  of  Nurick^  and  with  cavitation  frequency  measurements 
of  Chandra  and  Collicott^^.  Figure  3  provides  a  comparison  of  the  “critical  pressure  ratio”  measured  by 
Bergwerk  and  the  results  from  the  model.  Excellent  agreement  is  shown  over  the  range  of  L/jD  values 
presented. 

Comparisons  on  the  discharge  coefficient  (Fig.  4)  showed  greater  discrepancies.  The  discharge  coefficient 
is  defined  as  the  ratio  of  the  measured  (or  calculated)  flowrate  to  the  ideal  flowrate  one  gets  by  using  the 
Bernoulli  velocity.  This  parameter  measures  contraction  and  friction  losses  in  this  particular  flowfield.  Since 
the  flow  is  inherently  unsteady  for  most  cases  studied,  we  computed  a  time-averaged  value,  Cd,  which  was 
compared  to  Nurick’s  measurements  in  Fig,  4.  The  model  consistently  under  predicts  consistently  under 
predicts  discharge  coefficient  values  by  a  few  percent  for  the  fairly  long  [L/D  =  6)  orifice  studied.  While 
the  calculations  assume  a  sharp-edged  inlet,  the  amount  of  inlet  rounding/beveling  on  the  experimental 
hardware  is  not  known.  Note  that  a  small  degree  of  rounding  can  lead  to  substantial  changes  in  discharge^® 
(also  see  Fig.  16  and  associated  discussion). 

The  dimensional  cavitation  frequency,  /* ,  is  compared  with  the  measurements  of  Chandra  and  Collicott^^ 
in  Fig.  5.  The  computed  results  were  obtained  by  comparing  2-D  simulations  with  measurements  on  a 
high  aspect  ratio  slot.  For  the  numerical  results,  cavitation  frequency  was  obtained  by  performing  Fourier 
transform.  Results  in  Fig.  5  compare  the  primary  harmonic  (which  contains  the  bulk  of  the  oscillatory 
energy)  for  both  experiment  and  calculation.  Good  comparisons  are  noted,  thereby  validating  the  utility  of 
the  model  to  assess  frequencies  developed  in  these  flows. 

Results  of  Parametric  Studies 

Given  the  nondimensionalization  describe  above,  a  sharp-edged  orifice  geometry  is  completely  prescribed  from 
an  input  L/D  value.  From  Eq.  1,  Reynolds  and  cavitation  numbers  also  characterize  the  flow.  However,  the 
flow  is  not  uniquely  defined  by  these  two  parameters  since  K  depends  on  the  pressure  drop,  vapor  pressure, 
and  inlet  pressure.  We  chose  to  use  the  downstream  pressure,  P2  as  the  third  parameter  characterizing  a 
given  flow  situation.  By  raising  P2  and  maintaining  a  constant  pressure  drop,  pressure  levels  throughout  the 
orifice  are  increased  and  the  extent  of  cavitation  tends  to  be  reduced. 

In  this  work,  we  focus  on  the  unsteadiness  caused  by  instability  of  the  vena-contracta  under  cavitating 
and  non-cavitating  conditions.  Previous  works^®”^^  provide  in-depth  maps  of  the  flowfield;  the  present  study 
focuses  primarily  on  unsteadiness  in  global/integral  parameters.  Nearly  200  simulations  were  conducted  in 
the  studies.  A  typical  simulation  required  6-8  hours  runtime  on  a  450  MHz  Pentium  II  PC. 

Figure  6  presents  the  gross  structure  of  the  flowfield  near  the  inlet  lip  at  a  given  instant  in  time  for 
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both  cavitating  and  non-cavitating  flows.  The  radial  coordinate  is  expanded  in  these  plots  for  clarity.  The 
pressure  increases  at  the  aft  end  of  the  vena  contracta  are  fed  forward  and  effect  the  subsequent  shape  of 
this  region.  This  process  is  documented  in  detail  for  cavitating  flows  in  Ref.  21.  In  the  case  where  cavitation 
is  present,  the  recirculation  zone  can  detach  from  the  main  vena  contracta  and  be  convected  downstream. 
The  particular  image  in  Fig.  6  depicts  an  instant  in  this  process.  The  influences  on  orifice  discharge 
characteristics  are  investigated  with  respect  to  this  unsteadiness  from  these  vena  contracta  instabilities  in 
the  following  sections. 

Orifice  Flow  Under  Partial  Cavitation  Conditions 

A  “baseline”  case  is  presented  to  provide  the  reader  with  insight  into  the  unsteady  massflow  characterise 
tics  which  can  be  attributed  to  partial  cavitation.  Orifice  massflow  is  represented  in  terms  of  a  discharge 
coefficient,  Cdj  which  is  defined  in  the  usual  manner: 

CD=ni/{pivA)  (7) 

where  m  is  the  massflow  obtained  by  integration  of  the  velocity  profile  at  the  exit  plane,  v  is  the  Bernoulli 
velocity  as  described  previously,  and  A  is  the  orifice  cross-sectional  area.  The  conditions  noted  in  Figs.  7 
and  8  would  be  consistent  with  the  flow  of  water  through  a  580  micron  orifice  with  an  inlet  pressure  (Pi)  of 
4.3  atmospheres. 

Figure  7  shows  the  periodic  variations  in  dimensionless  cavitation  length  {Lc)  for  these  conditions.  In 
this  work,  the  cavitation  length  was  assumed  to  correspond  to  the  downstream  point  in  the  orifice  where  the 
pseudo-density  had  a  value  of  0.98.  Since  the  density  changes  rapidly  at  the  end  of  the  cavitation  region, 
similar  results  are  obtained  for  other  density  threshold  values.  Note  the  periodic  behavior;  in  this  case,  the 
cavity  oscillates  over  a  range  of  about  10-50%  of  the  orifice  length. 

The  repercussions  of  this  oscillatory  behavior  on  orifice  massflow  characteristics  are  shown  in  Fig.  8. 
In  this  particular  case,  the  unsteadiness  leads  to  orifice  massflow  variations  of  roughly  ±4  percent.  Most 
current  atomization  theories  suppose  that  the  flow  from  the  orifice  is  steady  and  that  the  fluid  exits  the  orifice 
initially  as  an  undisturbed  cylinder  subject  to  aerodynamic  interactions  from  the  gas  phase.  In  the  present 
case,  the  variations  in  massflow  would  lead  to  high- amplitude  perturbations  with  respect  to  a  linearized 
theory  such  as  this.  Experimentalists  generally  measure  massflow  by  collecting  fluid  over  some  time  period 
and  therefore  only  measure  an  averaged  discharge  coefficient.  For  this  reason,  the  unsteadiness  in  orifice  flow 
has  not  been  studied  in  any  detail  previously. 

By  performing  a  Fourier  transform  on  the  Lc  and  Cd  histories,  one  can  obtain  the  principle  frequency 
of  oscillation  of  each  of  these  parameters,  fc  and  fj^.  Since  these  frequencies  are  nondimensionalized  by  the 
orifice  diameter  and  the  Bernoulli  velocity,  they  represent  Strouhal  numbers  for  each  of  the  oscillations.  In 
addition,  we  should  point  out  that  the  Fourier  transforms  typically  show  most  of  the  energy  in  2-3  discrete 
frequencies^ The  bulk  of  the  energy  is  contained  in  the  primary  harmonic,  other  frequencies  which 
appear  to  be  higher  harmonics  tend  to  contain  much  less  energy.  This  work  will  focus  on  behavior  of  the 
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primary  harmonic  found  for  each  particular  orifice  flow. 

Figure  9  shows  the  variation  in  non-dimensional  fc  and  fo  with  variations  in  K  at  fixed  Reynolds 
number  for  an  orifice  with  L/D  S.  The  point  of  cavitation  inception  is  shown  on  the  figure;  the  region  to 
the  left  of  this  point  is  where  cavitation  is  evident  in  the  flowfield.  Over  this  range  note  the  agreement  in 
the  two  frequencies;  this  factor  confirms  the  notion  that  cavitation  oscillations  are  responsible  for  massflow 
variations  when  cavitation  is  present.  Note  that  cavitation  does  tend  to  increase  frequencies  up  to  a  point, 
beyond  which  the  very  large  cavities  tend  to  respond  at  lower  frequencies.  This  behavior  is  consistent  with 
prior  simulations^^  which  showed  the  oscillation  generated  by  a  disturbance  in  the  collapse  zone  which  is  fed 
upstream  at  a  roughly  constant  velocity. 

It  may  be  possible  to  make  use  of  this  result  to  design  passive  oscillations  of  a  desired  frequency.  In 
general,  the  frequencies  are  consistent  with  the  time  it  takes  a  fluid  element  to  traverse  the  length  of  the 
orifice  passage.  Atomizers  which  produce  fine  sprays  tend  to  produce  droplets  at  much  higher  frequencies 
than  those  noted  in  Fig.  9.  Partial  cavitation  may  be  used  to  “pump”  instabilities  for  atomizers  operating 
at  more  modest  pressure  drops  such  as  a  flow  produced  in  inkjet  printing  applications. 

Reynolds  Number  Effects 

A  series  of  simulations  were  conducted  to  investigate  the  influence  of  Reynolds  number  as  an  independent 
parameter.  Orifice  geometry,  cavitation  number,  and  exit'  pressure  were  all  held  fixed  in  this  series  of  runs. 
Results  on  the  mean  and  oscillatory  components  of  Lc  and  Cd  are  shown  in  Fig,  10.  The  mean  discharge 
coefficient  {Cd)  was  unaffected  by  changes  in  Re.  This  behavior  is  attributed  to  the  fact  that  a  relatively 
short  orifice,  L/D  =  4,  was  used  such  that  contraction  losses  dominated  the  losses  and  wall  friction  was  of 
lesser  importance.  As  Re  was  increased,  the  average  cavitation  length,  Lc,  also  changed  very  little.  However 
at  low  Re^  the  extent  of  cavitation  was  reduced;  below  Re  of  about  5000,  no  cavitation  was  observed  in  this 
case. 

The  RMS  oscillations  in  cavitation  length  and  discharge  coefficient  (ALc,  ACd)  show  a  similar  Insensi¬ 
tivity  to  Re  above  Re  ^  20, 000.  Below  Re  ^  5000  the  oscillations  vanish  and  a  steady  flow  is  computed.  It 
is  interesting  to  note  that  this  threshold  is  near  experimental  values  quoted  for  transition  to  turbulent  flow 
in  pipes.  In  this  problem,  the  unsteadiness  which  has  traditionally  been  attributed  to  turbulence  may  in 
many  cases  be  due  to  instability  of  the  vena-contract  a  since  our  laminar  calculations  show  “turbulent  like” 
qualities. 

Supply  Pressure  Effects 

A  series  of  simulations  were  conducted  with  a  fixed  orifice  diameter,  discharge  pressure,  and  fluid  properties 
to  identify  the  influence  of  supply  pressure  on  the  flowfield.  Figure  11  highlights  massflow  variations  for 
three  different  orifice  lengths  as  a  function  of  supply  pressure.  The  point  where  cavitation  begins  is  noted  in 
each  curve;  the  region  to  the  right  of  this  point  corresponds  to  partially  cavitated  flows.  Note  that  cavitation 
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becomes  evident  at  supply  pressures  in  the  2-3  atmosphere  range  for  these  sharp-edged  inlets;  since  most 
applications  utilize  higher  supply  pressures  cavitation  is  presumed  to  be  present  in  a  large  fraction  of  pressure 
atomizers  in  use  today.  The  highest  massflow  variations  occur  in  this  region  as  the  violence  associated  with 
cavity  collapse  (or  partial  collapse)  and  reformation  lead  to  significant  changes  in  massflow.  It  is  interesting 
to  note  that  significant  massflow  variations  are  shown  even  for  non-cavitating  flows.  Once  again,  these 
variations  are  largest  for  the  shortest  orifice  passages. 

Figure  12  shows  the  oscillations  in  cavitation  length  with  increased  supply  pressure.  The  increases  in 
ALc  are  directly  correlated  with  increases  in  ACd  as  discussed  previously.  Note  that  as  the  orifice  length 
increases,  ALc  becomes  relatively  insensitive  to  this  parameter,  the  cavitation  is  restricted  to  a  fairly  small 
fraction  of  the  orifice  length. 

Figure  13  depicts  the  changes  in  non-dimensional  oscillation  frequency  with  Pi  for  each  of  the  three 
orifices  studied.  There  is  a  general  trend  toward  decreasing  frequencies  with  increased  orifice  length;  the 
frequencies  are  coupled  to  the  orifice  transit  time.  The  onset  of  cavitation  signals  a  dramatic  increase  in 
frequency  for  each  case  studied.  It  may  be  possible  to  use  this  fact  as  a  diagnostic  tool  by  instrumenting 
an  orifice  with  microphones  or  transducers  to  detect  acoustic  energy  in  the  fluid.  In  combining  these  results 
with  Fig.  11  one  can  note  that  peak  oscillation  frequencies  are  not  correlated  with  maximum  variations  in 
massflow;  instead  massflow  variations  are  correlated  with  Lc^ 

Effect  of  Discharge  Pressure 

Figure  14  highlights  the  effect  of  discharge  pressure  (P2)  on  the  range  under  which  cavitation  is  prominent 
in  a,  L/D  —  6  orifice.  Increasing  back  pressure  leads  to  modest  increases  in  the  K  value  at  which  inception 
occurs.  The  flipping  condition  appears  to  be  insensitive  to  back  pressure,  occurring  near  K  =  1.5  for  the 
three  cases  investigated.  Of  course,  the  flipping  condition  is  highly  dependent  on  orifice  length. 

Finally,  we  note  that  the  cavitation  range  tends  to  be  extended  as  the  discharge  pressure  is  raised.  High 
pressure  applications  such  as  combustion  chambers  are  subject  to  cavitating  flows  over  a  much  wider  range  of 
pressures  than  one  might  see  in  an  ambient  pressure  test.  This  points  out  the  scaling  problems  encountered 
in  cold  flow  testing  of  high  pressure  combustion  systems;  cavitation  and  Reynolds  numbers  cannot  both  be 
matched  unless  both  pressure  levels  and  pressure  drops  are  replicated  in  testing. 

The  implications  of  changes  in  back  pressure  on  oscillation  frequencies  is  shown  in  Fig.  15.  Raising 
pressure  levels  in  the  atomizer  leads  to  increases  in  oscillation  frequencies  over  the  range  of  conditions  in 
which  cavitation  is  present.  This  behavior  is  attributed  to  the  fact  that  pressure  gradients  increase  at  constant 
K  with  the  increase  in  discharge  pressure  since  the  orifice  dimensions  are  fixed  in  this  particular  study.  Higher 
pressure  gradients  lead  to  higher  velocities  feeding  periodic  processes  associated  with  the  unstable  cavitation 
region.  Under  non-cavitating  conditions,  the  change  in  P2  had  no  effect  on  the  oscillation  frequencies. 
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Effect  of  Inlet  Rounding 

A  series  of  simulations  were  conducted  with  fixed  fluid  properties,  discharge  pressure,  and  L/D  io  assess  the 
effects  of  inlet  rounding.  Several  experiments  and  numerical  simulations^®  have  noted  that  inlet  rounding 
can  delay  the  occurrence  of  cavitation  and  increase  discharge  coefficient  through  reductions  in  contraction 
losses.  However,  the  effects  of  rounding  on  the  unsteady  behavior  has  not  been  investigated  in  any  great 
detail.  For  this  reason,  the  unsteady  characteristics  have  been  studied  assuming  the  inlet  is  rounded  with 
a  circle  of  given  radius.  Figure  16  depicts  the  level  of  massflow  variations  due  to  unsteadiness  at  various 
cavitation  numbers  for  three  levels  of  inlet  rounding  (r/ii  =  0,  0.05,  0.1).  The  cavitation  onset  point  is 
noted;  cavitating  results  lie  to  the  left  of  this  point  since  decreasing  K  corresponds  to  increased  supply 
pressure. 

Results  in  Fig.  16  show  that  massflow  variations  are  reduced  dramatically  for  even  a  small  amount  of  inlet 
rounding.  For  an  inlet  radius  of  10%  of  the  orifice  radius  {r/R  =  0.1)  massflow  variations  have  all  but  been 
eliminated  for  the  entire  pressure  range  noted.  These  results  are  disturbing  in  that  there  is  a  great  sensitivity 
of  the  flow  unsteadiness  to  minor  changes  in  shape.  Since  real  orifices  are  not  infinitely  sharp,  cavitation  sites 
will  most  likely  appear  at  burrs  or  notches  formed  during  the  fabrication  process.  These  features  are  generally 
not  reproducible  and  lead  to  hole-to-hole  variations  which  are  known  to  plague  some  manufacturers.  For 
example,  modern  diesel  orifices  are  typically  2-300  microns  in  diameter  and  defects/ tolerances  are  generally 
a  significant  fraction  of  this  diameter.  Reproducibility  in  tiny  orifices  can  be  problematic  for  these  reasons. 

Figure  17  highlights  cavitation  length  variations  from  the  inlet  rounding  study.  The  overall  level  of  length 
variations  is  reduced  as  inlet  rounding  is  increased.  In  addition,  rounding  of  the  inlet  suppresses  cavitation 
(and  ATc)  moving  inception  and  partial  cavitation  to  lower  K  and  higher  Pi  values.  It  is  interesting  to 
note  that  substantial  ALc  values  are  still  noted  for  the  rounded  inlets  (compare  values  with  those  in  Fig. 
12),  yet  the  repercussions  on  the  massflow  are  much  less  pronounced.  The  cavity  still  oscillates,  but  the 
undulations  are  much  less  violent  than  in  the  case  of  a  sharp-edged  inlet. 

Conclusions 

The  internal  flow  in  plain  orifice  atomizers  operating  at  modest  or  high  Reynolds  numbers  is  inherently 
unsteady  due  to  the  instability  of  the  vena-contract  a  formed  at  the  inlet  lip.  This  phenomenon  leads  to 
periodic  variations  in  orifice  massflow  over  a  time  scale  near  the  orifice  transit  time.  Cavitation  is  shown 
to  enhance  both  the  magnitude  and  the  frequency  of  these  oscillations.  The  influence  of  Reynolds  number, 
supply  pressure,  discharge  pressure,  orifice  L/D  and  inlet  rounding  have  been  investigated  using  laminar 
axisymmetric  calculations  with  a  homogeneous  flow  model  for  addressing  cavitation. 

Both  the  mean  and  oscillatory  components  of  cavitation  length  and  discharge  coefl6cient  are  affected 
little  by  Reynolds  numbers  when  it  is  above  about  20,000.  Steady  flow  solutions  were  found  for  Reynolds 
numbers  below  about  5000.  Unsteadiness  in  orifice  massflow  tended  to  increase  with  decreased  orifice  L/D. 
This  unsteadiness  also  increased  rapidly  near  conditions  where  the  orifice  was  cavitating.  Modest  amounts  of 


19 


cavitation  increased  the  frequency  of  massflow  oscillations  with  a  maximum  frequency  occurring  at  modest 
cavitation  lengths.  Increased  discharge  pressure  tended  to  delay  cavitation  inception  and  broaden  the  range  of 
injection  pressures  over  which  cavitation  was  present.  Rounding  of  the  inlet  had  dramatic  effects  in  reducing 
unsteadiness  in  massflow,  even  under  cavitated  conditions.  The  large  sensitivity  due  to  minor  geometry 
changes  is  unsettling  in  that  minor  manufacturing  variations  can  play  important  roles  in  the  overall  behavior 
of  these  atomizers. 


Nomenclature 

Cd  ~  Orifice  Discharge  Coefficient  (Eq.  7) 

D  -  Orifice  diameter 

/  -  Non-dimensional  frequency 

/*  -  dimensional  frequency 

K  -  Cavitation  number 

L  -  Orifice  length 

m  -  Orifice  massflow 

Lc  -  Cavitation  length 

n  -  Bubble  number  density 

P  -  Pressure 

Re  -  Reynolds  number 

r  -  Radial  or  transverse  coordinate 

t  -  Time 

V  -  Bernoulli  velocity 
z  -  Axial  coordinate 
a  -  Void  fraction 
fi  -  Viscosity 
p  -  Fluid  pseudo-density 


Subscripts 

1  -  Inlet 

2  -  Outlet 

C  -  Cavitation  length 
D  -  Discharge  coefficient 
/  -  Liquid 
V  -  Vapor 
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Figure  Captions 

1.  Schematic  of  the  flow  features  in  a  cavitating  orifice. 

2.  Computational  domain  and  mesh  for  orifice  flow  (for  clarity,  only  a  small  portion  of  the  grid  lines  are 
shown) . 

3.  Critical  pressure  ratio  comparison  with  experimental  results^;  D=2.5mm,  Re  ^  10000,  sharp-edged 
inlet. 

4.  Discharge  coefficient  Co  comparison  with  experimental  results®;  L/D=6,  D==3.18  mm,  P2~13.8  psi. 

5.  Cavitation  frequency  comparison  with  experimental  results^  ^  for  high  aspect  ratio  slots 

6.  Typical  streamline  patterns  near  the  inlet  lip  under  cavitating  {Re  =  16188,  A"  =  1.53)  and  noncavi- 
tating  {Re  =  11849,  K  =  1.98)  conditions.  The  vertical  scale  has  been  expanded  for  clarity. 

7.  Typical  quasi- periodic  behavior  of  cavitation  length 

8.  Typical  orifice  discharge  coefficient  {Cd)  behavior. 

9.  Frequency  of  discharge  (/o)  and  and  cavitation  length  {fc)  vs.  cavitation  number;  L/D=8,  P2  =  latm. 

10.  The  effect  of  Reynolds  number  on  average  and  unsteady  components  of  discharge  coefficient  and  cavi¬ 
tation  length;  L/D=4,  K  =  1.6,  P2=l  atm. 

11.  The  effect  of  supply  pressure  on  the  fluctuation  of  discharge  coefficient,  0.566  mm,  P2=l  atm, 
sharp-edged  inlet. 

12.  The  effect  of  supply  pressure  on  fluctuations  in  cavitation  length,  D— 0.566mm,  P2=l  atm,  sharp-edged 
inlet. 

13.  The  effect  of  supply  pressure  on  the  frequency  of  discharge  coefficient,  D— 0.566  mm,  P2“l  atm,  sharp- 
edged  inlet. 

14.  The  effect  of  back  pressure  on  the  flow  state  L/D=6,  D=0.566  mm. 

15.  The  effect  of  supply  pressure  on  the  frequency  of  discharge  coefficient,  D=0.566  mm,  P2— 1  atm,  sharp- 
edged  inlet. 

16.  The  effect  of  inlet  rounding  on  fluctuations  of  orifice  massflow,  L/D=6,  P2=l  atm. 

17.  The  effect  of  inlet  rounding  on  fluctuations  in  cavitation  length,  L/D=6,  P2=l  atm. 
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Figure  2:  Computational  domain  and  mesh  for  orifice  flow  (for  clarity,  only  a  small  portion  of  the  grid  lines 
are  shown). 


Figure  3:  Critical  pressure  ratio  comparison  with  experimental  results^;  D=2.5mm,  Re  ^  10000,  sharp-edged 
inlet. 
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Figure  4:  Discharge  coefficient  Cd  comparison  with  experimental  results^;  L/D:=6,  D=3.18  mm,  F2=13.8 


Figure  5:  Cavitation  frequency  comparison  with  experimental  results^ ^  for  high  aspect  ratio  slots 


Figure  6:  Typical  streamline  patterns  near  the  inlet  lip  under  cavitating  {Re  =  16188jA"  =  1.53) 
noncavitating  {Re  =  11849,  K  —  1.98)  conditions.  The  vertical  scale  has  been  expanded  for  clarity. 
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Figure  7:  Typical  quasi-periodic  behavior  of  cavitation  length  Lc 


Figure  8:  Typical  orifice  discharge  coefficient  (Co)  behavior. 
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Figure  10:  The  effect  of  Reynolds  number  on  average  and  unsteady  components  of  discharge  coefficient  and 
cavitation  length;  L/D=:4,  K  =  1.6,  P2=l  atm. 
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Figure  11:  The  effect  of  supply  pressure  on  the  fluctuation  of  discharge  coefficient,  D=0.566  mm,  atm, 
sharp-edged  inlet. 


Figure  12:  The  effect  of  supply  pressure  on  fluctuations  in  cavitation  length,  D=::0. 566mm,  P2^l  atm, 


Figure  16:  The  effect  of  inlet  rounding  on  fluctuations  of  orifice  massflow,  L/D=:6,  P2=l  atm. 


8  Appendix  B  -  2D  and  Axisymmetric  Turbulenct,  Cavitating 
Flow  Simulations 


Xu,  C.,  “Simulation  of  Orifice  Internal  Flows  Including  Cavitation  and  Turbulence”  Ph. 
D.  Thesis,  School  of  Aeronautics  and  Astronautics,  Purdue  University,  West  Lafayette,  IN, 
August  2001,  pp.  44  -  81.  1996. 
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4.  2D  AND  AXISYMMETRIC  TURBULENT,  CAVITATING  FLOW 

SIMULATIONS 


4.1  Introduction 

The  presence  of  cavitation  inside  a  slot /orifice  injector  nozzle  increases  the  tur¬ 
bulence  level  of  the  flow  greatly.  Recent  experiments  (Tamaki  et  ah,  1998;  Hiroyasu, 
2000)  have  shown  that  the  turbulence  in  the  nozzle  hole  resulting  from  cavitation 
is  a  mechanism  that  promotes  atomization.  Under  partial  cavitation  conditions,  the 
cavitation  region  extends  only  a  fraction  of  the  orifice  length,  its  shape  is  unsteady, 
and  oscillates  in  a  quasi-periodic  manner  (Bunnell  and  Heister,  1999;  Henry,  1997; 
Sanchez,  1999).  At  the  rear  end  of  the  cavity,  the  main  flow  of  the  liquid  eventually 
reattaches  to  the  wall,  and  clusters  of  bubbles  collapse  accompanying  the  recovery 
of  pressure.  This  process  generates  a  great  amount  of  turbulence  in  the  flow.  Fig. 
4.1  is  a  schematic  representation  of  this  process.  While  the  mechanism  of  energy 
transfer  from  mean  flow  to  turbulence  during  cavitation  remains  unresolved,  experi¬ 
ments  have  shown  that  turbulence  intensity  is  increased  in  the  presence  of  cavitation. 
Gopalan  and  Katz  (2000)  observed  that  the  unsteady  cavity  collapse  involves  sub¬ 
stantial  increases  in  turbulence  intensity,  momentum  and  displacement  thicknesses  in 
the  boundary  layer.  They  also  showed  that  the  collapse  of  bubbles  is  the  dominant 
source  of  vorticity  downstream  of  a  cavity.  Ruiz  and  He  (1999)  have  shown  that  there 
is  a  jump  in  turbulence  intensity  in  a  flow  after  a  cavitation  zone,  and  the  turbulence 
decays  more  slowly  than  ordinary  turbulence. 

In  this  chapter,  the  turbulent  characteristics  of  cavitating  flows  are  investigated. 
The  k  —  u  turbulence  model  and  Chen  and  Heister’s  (1995)  homogeneous  fluid  model 
are  employed  in  simulations  of  turbulent  cavitating  flows.  The  predicted  discharge 
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After  reattachment  accompanion 

with  bubble  detachment  and  collaoses.  Turbulence  enhances 


Figure  4.1  Partial  cavitation  condition  enhances  atomization 


coefficient,  cavitation  length  and  turbulence  quantities  are  compared  with  available 
experimental  data. 

4.2  Model  Formulation 

In  order  to  avoid  the  presence  of  additional  correlations  due  to  the  presence  of 
density  fluctuation,  it  is  convenient  to  recast  the  instantaneous  flow  fields  in  terms 
of  Favre  averages  instead  of  time  averages.  A  Favre  averaged  variable,  denoted  by  a 
tilde,  is  defined  as  _ 

/  =  ^  (4.1) 

P 

Here,  we  introduce  Favre  averaging  because  the  density  changes  greatly  between  the 
inside  and  outside  of  the  cavity  region.  Favre  averaging  is  a  mathematical  simpli¬ 
fication  to  eliminate  the  additional  terms  resulting  from  time-averaging  because  of 
the  fluctuation  of  density  (Wilcox,  1998).  Thus,  we  decompose  the  flow  properties  as 
follows 

Uj  —  Uj  -j- 

p  =  P  +p'  (4.2) 

P  =  P  +  p' 

where  is  Favre  mass-averaged  velocity,  and  u"  is  the  fluctuation  part  of  the  velocity. 
Density  and  pressure  are  decomposed  as  the  time-averaged  part  and  fluctuation  part. 
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The  governing  equations  are  as  shown  below. 


a(r»  d{r’^pui)  ^  ^ 


dt 


d{r^pui)  d{r‘^pujUi)  ^ 


dt 


dxi 


9P  a  9  rr 
T  ^  O  T 


dx 


dxi 


jii 


(4.3) 


(4.4) 


where,  a  =  0  for  a  2-D  problem  and  a  =  1  for  an  axisymmetric  problem,  tij  and  r,j 


are  the  viscous  stress  tensor  and  mass-averaged  Reynolds  stress  tensor  Tij  =  —pu'-u'-, 
respectively,  which  are  defined  as 


t  -  —  —  (  S- 
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Here,  the  viscous  and  Reynolds  stresses  are  expressed  in  compressible  form  since  the 
divergence  of  velocity  is  not  zero  in  the  cavity  region.  The  eddy  viscosity,  pj,  is 
calculated  from  the  k  —  to  model,  and  Sij  is  strain  rate  tensor. 


1  /  dui  duj 

2  dxi 


The  k  —  u)  model  (Wilcox,  1998)  takes  the  following  form 

du 
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(4.9) 


One  critical  issue  is  how  to  resolve  the  eddy  viscosity  in  the  two  phase  mixture  region. 
Here,  we  took  the  same  treatment  for  the  eddy  viscosity,  pj,  as  that  for  the  dynamic 
viscosity,  p,  for  two  phase  mixture, 

PT  oc  (l-fg)-  =p- 
CO  CO 


(4.10) 


37 


where  fg  is  void  fraction,  p  —  1  —  fg.  The  implication  of  Eq.  4.10  is  that  the  contribu¬ 
tion  of  gas  phase  to  turbulence  is  neglected  compared  to  that  of  liquid  phase  because 
pg  «  pi-  Including  the  closure  coefficient,  eddy  viscosity  is  evaluated  by 


Pt 


* — 
a  p- 


iO 


(4.11) 


The  closure  coefficients  and  auxiliary  relations  of  the  k  —  oj  model  are  given  in  sec¬ 
tion  3.2.  All  quantities  in  the  above  equations  are  non-nondimensionalized  based  on 
Bernoulli  velocity  U*,  slot  height /orifice  diameter  D*,  liquid  density  p]  and  dynamic 
viscosity  Hence 


u*  =  U*ui  P*  =  p]U*^P  X*  =  D*x 


k*  =  U*^k 


U*2 


-U 


l^l 


Re  — 


PlU*D* 

Pi 


(4.12) 


4.3  Injector  Internal  Flows 

4.3.1  Boundary  Conditions 


For  flow  through  an  injector  slot/orifice,  a  constant  pressure  condition  is  utilized 
both  at  the  inflow  and  outflow  boundaries.  To  insure  the  accuracy  of  the  constant 
pressure  inflow  condition,  the  inflow  boundary  is  placed  five  gap  heights  or  orifice 
diameters  upstream  of  the  inlet  corner.  Previous  simulations  (Bunnell  et  al.,  1999) 
used  zero  velocity  at  the  inflow  boundary.  Because  the  inflow  boundary  is  far  from 
the  inlet,  and  the  velocity  there  approaches  zero,  this  only  leads  to  minor  errors  in 
the  calculation.  However,  this  treatment  violates  conservation  of  mass  flow,  since 
zero  velocity  at  the  inflow  means  zero  mass  flow  rate  at  the  inflow.  In  order  to 
circumvent  this  dilemma,  we  employ  a  sink  at  the  origin  to  approximate  the  inflow 
velocity  conditions.  A  similar  treatment  was  utilized  before  to  approximate  the  inflow 
conditions  of  a  three  dimensional  manifold  cross-flow  by  Bunnell  (1999).  As  shown 
in  Fig.  4.2  an  artificial  sink  is  put  at  the  origin.  For  a  2-D  problem  the  velocity  at 
the  inflow  boundary  is  calculated  by 
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Figure  4.2  Schematic  representation  of  a  sink  approximation  of  inflow  velocity 

The  strength  of  the  sink,  A,  is  updated  during  each  time  step  by  the  conservation 
of  mass  flow  rate  through  the  nozzle  passage.  The  only  diff'erence  between  an  ax- 
isymmetric  problem  and  a  2-D  problem  at  this  point  is  that  a  three  dimensional  sink 
is  utilized  for  the  former.  This  assumption  might  be  not  fully  proper  for  unsteady 
flow  conditions.  However,  the  unsteadiness  is  considerably  smaller  than  the  average 
magnitude  of  mass  flow  rate  through  the  nozzle  hole.  This  treatment  also  results  in 
the  improvement  of  discharge  coefficient  predictions  which  will  be  discussed  in  section 
4.3.3. 

The  inflow  boundary  condition  of  turbulence  properties  k  and  uj  have  an  un¬ 
expected  effect  on  the  turbulence  models’  performance.  The  profile  for  k  can  be 
estimated  from  experimental  data.  However,  it  is  hard  to  set  the  inflow  condition  for 
u.  Even  if  k  and  are  sufficiently  small  at  the  inflow  boundary,  the  choice  of  ui  can 
have  a  significant  effect  on  the  flow  downstream.  It  has  been  documented  that  free 
shear  flow  spreading  rates  are  sensitive  to  the  free  stream  value  of  u>  (Wilcox,  1998). 
In  the  flow  over  a  backward  facing  step,  the  prediction  gives  a  shorter  reattachment 
length  with  a  larger  initial  length  scale  (larger  eddy  viscosity)  (Nallasamy,  1987).  For 
this  particular  problem,  the  inflow  boundary  is  put  five  gap  heights  or  orifice  diame¬ 
ters  upstream  of  the  inlet  corner  where  the  velocity  is  small.  Because  no  experimental 
data  are  available  and  there  is  no  reason  to  set  a  high  k  value  there,  we  choose  a  k 
value  of  0.5%t/j^  at  the  inflow  boundary,  which  is  very  small  due  to  the  small  inlet 
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Laminar  calculation  Turbulent  calculation 


Figure  4.3  Sensitivity  of  cavitation  length  to  inflow  uj  conditions 


velocity.  The  specific  dissipation  rate  at  the  inflow  boundary  is  given  by 

^/k 


(4.14) 


A  is  a  scale  factor  and  can  be  a  free  parameter  to  adjust  the  inflow  length  scale.  We 
performed  a  study  of  the  effect  of  inflow  u;  value  on  the  flow  downstream.  In  Fig.  4.3, 
the  curve  on  the  left  is  from  a  laminar  calculation  which  represents  zero  length  scale 
at  the  inflow,  the  curves  on  the  right  are  due  to  turbulent  calculations  with  three 
different  length  scales.  The  laminar  solution  gives  a  cavitation  length  with  significant 
oscillations  representing  strong  unsteadiness  in  the  flow  field.  However,  turbulent 
modeling  produces  a  constant  value  of  cavitation  length.  This  is  not  outside  our 
expectation  since  Reynolds  averaged  models  decompose  the  instantaneous  flow  field 
into  an  average  part  and  a  fluctuating  part  and  only  solves  for  the  mean  flow  part. 
This  simplification  usually  results  in  a  steady  state  solution  which  still  provides  the 
essence  of  a  problem.  In  this  study,  we  choose  three  different  A  values  75,  150  and 
300.  Larger  A  values  correspond  to  smaller  a;  values  (larger  eddy  viscosity).  As  shown 
in  Fig.  4.3,  there  is  no  difference  in  the  resulting  Lc  when  A  =  150, 300.  Note  that 
the  cavitation  length  in  these  two  cases  is  close  to  the  average  value  of  Lc  from  the 
laminar  result.  We  prefer  a  smaller  A  value  which  gives  a  smaller  eddy  viscosity  at 
the  inflow  because  we  believe  the  turbulence  level  is  not  high  there.  A  A  value  of  150 
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Figure  4.4  Velocity  profile  resolved  by  different  meshes  (a) 

gives  an  inflow  eddy  viscosity  approaching  unity,  and  this  value  is  used  through  all 
the  calculations. 

4.3.2  Grid  Refinement  Study 

In  previous  work  (Bunnell  et  ah,  1999),  a  grid-refinement  study  was  performed. 
A  grid  of  140x60  was  found  to  be  fine  enough  for  laminar  simulations.  However  this 
grid  is  not  proper  for  turbulence  modeling  calculations.  It  is  necessary  to  place  a 
very  fine  grid  in  the  near  wall  region  to  resolve  a  turbulent  boundary  layer.  Since 
the  grid  used  is  highly  stretched  in  the  wall  normal  direction,  we  carried  out  this  grid 
dependency  study  in  two  ways:  (1)  the  effect  of  minimum  distance  of  the  grid  beyond 
the  wall,  (2)  the  effect  of  the  total  number  of  grid  points  in  the  normal  direction. 
Fig.  4.4  shows  the  effect  of  the  minimum  distance  away  from  the  wall  of  a  mesh  on 
the  velocity  profile  at  the  exit  of  a  circular  orifice.  Here  and  are  defined  as 

u'^  =  —  and  y'^  = -  (4-15) 

Ur  u 

where  Ur  is  the  friction  velocity. 

A  typical  velocity  profile  for  a  turbulent  boundary  layer  consists  of  three  regions: 
viscous  layer,  log  layer  and  defect  layer.  It  is  clearly  shown  in  Fig.  4.4  that  a  mesh 
with  the  first  point  placed  at  0.002  can  not  resolve  the  three-layer  structure.  As  the 
first  grid  moves  closer  to  the  wall  the  friction  velocity  Ur  resolved  is  greater,  generating 
a  lower  curve  in  the  plot.  However  a  mesh  with  Aymin  =  0.000125  produces  a  lower 
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Figure  4.5  Velocity  profile  resolved  by  different  meshes  (b) 

curve  than  the  one  produced  by  a  mesh  with  ^y-min  —  0.000080.  This  is  because 
a  constant  total  number  of  normal  grid  points  of  60  was  used  in  the  calculations  of 
Fig.  4.4.  The  smaller  Aymin  is,  the  coarser  the  grid  in  the  middle  of  the  flow  field.  To 
a  certain  degree,  such  a  grid,  however,  is  less  sufficient.  Fig.  4.5  shows  the  effect  of  the 
total  number  of  grid  points  in  the  normal  direction.  With  Aymin  =  0.000125  a  mesh 
of  140x90  is  equally  sufficient  compared  to  a  mesh  of  140x120  and  a  mesh  of  140x90 
with  Aymin  =  0.000080.  Thus  we  consider  a  grid  of  140x90  with  Aymin  —  0.000125 
to  be  fine  enough  for  the  turbulence  calculation.  With  such  a  choice,  the  wall  unit  of 
the  first  point  below  or  above  the  wall  is  at  j/+  =  0.16  and  about  10  points  are  placed 
below  =  2.5,  as  recommended  by  Wilcox  (1998). 

4.3.3  Discharge  Coefficient  Prediction 

A  few  axisymmetric  runs  were  made  to  compare  with  the  discharge  coefficient, 
Co,  measurements  by  Nurick  (1976)  on  a  circular  orifice.  Fig.  4.6  shows  the  results 
of  a  comparison  on  an  orifice  of  L/D  =  6,  D  =  3.18mm  under  a  back  pressure 
of  P2  =  13.8psi.  In  Fig.  4.6,  “Laml”  represents  a  laminar  calculation  with  zero 
inflow  velocity  boundary  condition,  “Lam2”  represents  a  laminar  calculation  with 
an  approximated  inflow  velocity  from  an  artificial  sink  as  described  in  the  previous 
section,  and  the  line  denoted  as  “Turb”  is  calculated  from  the  turbulence  model.  As 
shown  in  the  plot,  the  new  inflow  velocity  treatment  greatly  improves  the  prediction 


42 


K 


Figure  4.6  Discharge  coefficient  Cd  comparison  with  experimental  results; 

LJ D  =  Qj  D  =  3.18mm 


Figure  4.7  Comparison  of  velocity  profile  at  exit  for  laminar  and  turbulent 
solutions;  L( D  =  Q,  D  =  3.18mm 
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of  discharge  coefficient.  There  is  an  increase  in  the  magnitude  of  Cd  of  4.7%  over  the 
zero  inflow  velocity  treatment.  The  turbulence  model  shows  a  further  improvement 
and  gives  results  somewhat  closer  to  the  experimental  data.  The  differences  between 
the  turbulence  model  and  laminar  predictions  on  Cd  might  be  explained  by  Fig.  4.7 
in  which  the  exit  velocity  profiles  are  plotted  for  both  calculations.  The  turbulent 
velocity  profile  is  fuller  than  the  laminar  one.  Due  to  this  feature,  the  turbulence 
model  yields  a  larger  Cd  by  1.6%.  In  the  calculations  a  sharp-edged  inlet  is  assumed, 
however  the  amount  of  inlet  rounding  on  the  experimental  hardware  is  not  known.  As 
noted  in  section  2.4.5  small  degree  of  inlet  rounding  can  lead  to  substantial  increase 
in  discharge  coefficient.  This  factor  might  be  the  primary  reason  for  the  model’s 
consitent  underprediction  of  discharge  coefficient. 

4.3.4  Cavitation  Extent  Comparison 

The  2-D  code  was  used  to  simulate  the  flow  through  slots  with  the  geometry  used 
by  Henry  (1997)  and  Sanchez  (1999).  Extensive  comparisons  between  cavitation 
length  predicted  by  the  laminar  code  (Bunnell  et  al.  1999)  and  the  measurements 
provided  by  Henry  have  been  conducted.  The  measured  cavitation  length  falls  within 
the  range  of  the  maximum  and  minimum  cavitation  lengths  predicted  by  the  laminar 
code  for  various  slot  size  and  flow  conditions.  To  validate  the  turbulent  code,  a  similar 
comparison  between  the  average  cavitation  length  predicted  by  the  turbulent,  code 
and  the  measurements  by  Henry  and  Sanchez  was  performed.  The  computational 
conditions  are  presented  in  Table  4.1. 

As  mentioned  in  the  previous  section,  the  turbulent  code  produces  a  constant 
value  of  the  cavitation  length.  Figures  4.8-  4.10  show  the  results  for  the  three  slots. 
The  best  agreement  was  obtained  with  the  slot  with  medium  length/height  ratio, 
L/D=9.375.  For  this  case  the  predicted  average  cavitation  length,  Lc,  agrees  very  well 
with  Henry’s  data.  For  all  three  kinds  of  slots,  Henry  had  a  larger  cavitation  region 
than  Scanchez  had,  and  the  slot  with  the  largest  length/height  ratio,  L/D  =  10.714 
shows  the  largest  difference.  The  reason  for  the  variation  between  Henry  and  Sanchez 
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Table  4.1  Computational  conditions 


case 

AP(psi) 

Re 

K 

Pv 

L/D= 

=10.174 

D(mm): 

=0.889 

C(mm)=0.325 

1 

30 

22420 

1.48 

-0.239 

2 

32 

22420 

1.45 

-0.224 

3 

34 

23110 

1.42 

-0.211 

4 

36 

23780 

1.40 

-0.199 

5 

38 

24431 

1.38 

-0.189 

L/D: 

=9.375 

D(mm): 

=1.016 

C(mm)=0.325 

6 

28 

23968 

1.51 

-0.256 

7 

30 

24809 

1.48 

-0.239 

8 

32 

25623 

1.45 

-0.224 

9 

34 

26411 

1.42 

-0.211 

10 

36 

27177 

1.40 

-0.199 

11 

38 

27922 

1.38 

-0.189 

L/E 

1=6.25 

D(mm): 

=1.524 

C  (mm)  =0.325 

12 

22 

31868 

1.65 

-0.326 

13 

24 

33285 

1.60 

-0.299 

14 

26 

34644 

1.55 

-0.276 

15 

28 

35952 

1.51 

-0.256 

16 

30 

37214 

1.48 

-0.239 

L 

V _ 1 

45 


AP(psi) 


Figure  4.8  Cavitation  length,  Lc  comparison  with  experimental  results, 

LjD  =  10.714 


AP(psi) 


Figure  4.9  Cavitation  length,  Lc  comparison  with  experimental  results, 

LjD  =  9.375 
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Figure  4.10  Cavitation  length,  Lq  comparison  with  experimental  results, 

LjD  =  6.25 


Figure  4.11  Overlays  of  numerical  results  and  photographic  time  shot  of  cavitation 
region,  (AF*  ~  3Qpsi,  Lj D  =  10.714).  Here,  p  =  0.98  on  the  outermost  contour  with 

gradients  of  0.2  between  contours 
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data  is  not  clear;  however,  they  used  different  methods  for  measuring  the  cavitation 
length.  The  computational  Lp  for  LfD  —  10.714  slot  lies  in  between  Henry’s  results 
and  Sanchez’s  results  and  is  closer  to  Henry’s.  The  Lj D  =  6.25  slot  shows  the  poorest 
agreement  between  computed  results  and  experimental  measurements.  The  reason 
for  that  might  be  that  the  collapse  and  reformation  of  the  bubble  cavity  tends  to 
be  more  violent  for  a  shorter  slot,  and  the  turbulence  model  does  not  account  for 
the  turbulence  caused  by  bubble  collapse.  Although,  the  turbulence  model  smears 
out  the  oscillations,  it  generates  a  cavitation  length  in  acceptable  agreement  with  the 
experimental  data. 

The  turbulent  calculation  also  produces  a  single  contiguous  cavitation  region  near 
the  inlet  corner,  consistent  with  experimental  observation.  Fig.  4.11  shows  density 
contours,  which  denote  the  cavity  region,  obtained  from  laminar  (the  upper  one) 
and  turbulent  (the  lower  one)  calculations,  overlaying  one  photographic  snapshot  by 
Henry.  Although  the  laminar  simulation  results  in  an  overall  cavitation  extent  consis¬ 
tent  with  experiment,  it  indicates  two  separate  regions  of  cavitation.  The  turbulence 
model  improves  on  this  point  by  generating  a  single  cavitation  region  which  appears 
to  be  quite  consistent  with  experimental  results  both  in  axial  and  cross-stream  direc¬ 
tions. 

4.3.5  Velocity  Fields  and  Boundary  Layer  Thicknesses 

Fig.  4.12  shows  the  velocity  vectors  of  the  three  cases  [Lj D  =  10.714)  described 
in  Table  4.1.  The  broken  lines  in  this  plot  are  density  contours  of  a  value  0.98.  All 
three  cases  indicate  that  the  flow  reattaches  to  the  wall  at  the  rear  end  of  cavitation 
region.  The  close-up  of  the  streamlines  near  the  inlet  corner  of  the  case  AP  =  32psi 
is  shown  in  Fig.  4.13.  As  seen  from  this  figure  the  flow  separats  at  the  second  corner 
of  the  bevel. 

Fig.  4.14  shows  velocity  profiles  inside  the  location  of  the  cavity  region  for  cavitat- 
ing  and  non-cavitating  conditions.  The  flow  inside  the  cavitation  region  in  general  is 
slower  under  cavitating  conditions  than  under  non-cavitating  conditions.  Especially 
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Figure  4.14  Velocity  profiles  under  cavitating  and  non-cavitating  conditions, 
Re  =  23780,  ir  =  1.4,  L/D  =  10.714 
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Figure  4.15  Comparison  of  boundary-layer  thickness  5,  displacement  thickness  S*, 
momentum  thickness  0,  and  form  factor  H,  distribution  on  the  upper  wall  of  a 
nozzle  passage,  LjD  =  6.25,  Re  =  34644,  K  =  1.55 
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in  the  middle  of  the  cavity  (a;  =  2),  a  strong  reverse  flow  occurs  between  the  wall  and 
cavity.  At  the  end  of  the  cavitation  region  {x  =  3)  the  velocity  profiles  approach  the 
same  shape  for  both  cavitating  and  non-cavitating  conditions. 

As  shown  in  Fig.  4.15  cavitation  has  a  significant  effect  on  the  downstream  bound¬ 
ary  layer  thickness,  S,  displacement  thickness  5*,  momentum  thickness  6.  They  are 
defined  as 


D  y|u=0.99ur, 
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(4.16) 


D~8 

D 


«  =  E 


u 


D-S 


Uri 


1  - 


u 


Ur 


Ay 


where  Umax  is  the  maximum  u  value  in  a  specific  cross  section.  The  cavitation  region 
corresponding  to  the  situation  of  Fig.  4.15  extends  to  x  =  1.12.  Fig,  4.15  shows  the 
distribution  of  5,  5*,  6,  and  shape  factor  H  along  the  streamwise  direction  for  both 
cavitating  and  non-cavitating  conditions.  Even  though  the  turbulence  model  predicts 
S  distributions  which  are  almost  identical  for  both  cavitating  and  non-cavitating  con¬ 
ditions,  it  produces  an  increased  displacement  and  momentum  thicknesses  for  cav¬ 
itating  conditions.  Previous  researchers  (Gopalan  and  Katz,  2000;  Kubota  et.  al, 
1992)  have  observed/predicted  the  same  behavior. 

Fig.  4.16  shows  the  contour  of  vorticity  which  is  defined  as 


dv  du 


0Jz  =  - 


(4.17) 


dx  dy 

Note  the  maximum  and  minimum  value  of  uj  occurs  near  the  walls.  To  highlight  the 
interior  distribution  of  u  the  contour  levels  have  been  reduced  to  such  as  indicated 
in  the  plot.  As  seen  from  the  plot,  high  magnitude  vorticity  exist  downstream  of  the 
inlet  corner  where  strong  recirculation  occurs.  The  plot  also  indicates  a  big  portion  of 
this  internal  flow  remains  irrotational  with  vorticity  close  to  zero.  Fig.  4.17  shows  a 
detail  comparison  of  the  distribution  of  vorticity  in  normal  direction  at  different  axial 
locations.  As  the  flow  moves  toward  the  exit,  the  maximum  value  of  the  vorticity  in 
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Figure  4.16  The  distribution  of  vorticity,  LjD  =  6.25,  Re  —  34644,  K  =  1.55 


the  upper  wall  boundary  layer  reduces  greatly,  however  considerable  level  of  vorticity 
remains  at  the  exit.  Note  the  vorticity  level  at  the  edge  of  the  boundary  layer  of  the 
upper  wall  is  about  0.1. 


4.3.6  TKE  and  Reynolds  stresses 


We  present  the  contour  plots  of  TKE  (turbulent  kinetic  energy)  in  figures  4.18, 
4.19  and  4.20  for  three  operating  conditions  to  demonstrate  the  effect  of  cavitation 
on  turbulence.  The  corresponding  pseudo  density  and  eddy  viscosity  contours  for 
each  case  are  also  shown  in  the  figures.  As  seen  from  the  plots,  the  peak  region  of 
TKE  locates  right  at  the  rear  end  of  cavity  region,  the  wake  of  the  cavity,  for  all 
operating  conditions.  It  is  well  known  that  the  wake  where  bubble  detachment  and 
collapses  accompanied  with  pressure  recovery  occurs  is  the  most  ’turbulent’  region.  It 
is  interesting  to  note  that  the  calculation  can  capture  this  important  phenomenon  in 
cavitating  flows.  Also  note  that  as  bubble  cavity  becomes  bigger  (cavitation  number 
getting  lower)  the  location  of  TKE  peak  shifts  downstream  and  the  peak  value  of 
TKE  is  nearly  doubled  for  the  modestly  cavitating  case  (Fig.  4.20)  compaxed  with 
the  slightly  cavitating  case  (Fig.  4.18).  The  distribution  of  eddy  viscosity  shows 
similar  behavior  to  that  of  TKE.  Eig.  4.21  shows  the  profiles  of  TKE  at  the  exit  for 

However, 


h  Ul  APi 


these  three  cases.  One  might  expect  k  oc  (x  AP  or  7-^  =  =  *  r.  • 

®  ^  k2  Ul  AP2 

the  increase  in  k  is  larger  than  what  is  expected  based  on  the  increased  pressure  loss 
alone.  The  ratio  of  the  pressure  loss  AP  =  SQpsi  to  AP  =  2Spsi  is  1.29,  whereas 
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the  ratio  of  the  maximum  TKE  at  the  exit  of  these  two  cases  is  1.71.  This  indicates 
cavitation  plays  a  role  in  the  increase  of  k  value. 

The  2-D  code  was  also  used  to  analyze  the  flow  through  a  slot  with  the  geometry 
consistent  with  the  one  used  by  Ruiz  and  He  (1999).  In  this  case  a  sharp  cornered, 
T/D  =  6  and  D  =  25mm,  slot  was  used.  Fig.  4.22  shows  turbulence  intensities  in  the 
streamwise  and  normal  direction  and  Reynolds  shear  stress  compared  with  experimen¬ 
tal  data  (Ruiz  and  He,  1999)  at  different  streamwise  locations.  Here,  Rxx  =  yj pu"u" 
and  Ryy  =  yj pv"v",  representing  turbulence  intensity  in  x  and  y  directions  are  the  di¬ 
mensionless  mass  weighted  velocity  fluctuation  in  streamwise  and  normal  directions  , 
and  Rxy  —  pu"v"  is  the  mass  averaged  Reynolds  shear  stress.  The  cavitation  extent  is 
also  shown  in  this  plot.  In  the  experiment,  an  artificial  air  filled  cavity  was  introduced 
near  the  inlet  corner,  which  is  different  from  a  vapor  filled  cavity  generated  from  cav¬ 
itation.  In  the  calculation  the  outflow  pressure,  P2  was  used  as  a  free  parameter  (the 
lower  P2  is,  the  more  likely  the  flow  is  to  cavitate)  to  match  the  cavity  extent,  which 
is  about  one  fourth  of  the  slot  length.  The  Reynolds  number.  Re  =11280,  is  the  same 
as  that  in  the  experiment. 

The  turbulence  model  predicts  higher  values  of  turbulence  intensity  and  Reynolds 
shear  stress  for  cavitating  flow  conditions  than  for  non-cavitating  conditions.  Quan¬ 
titatively  the  calculation  results  are  comparable  to  the  experimental  data.  Near  the 

I  I 

rear  end  of  cavity  region  the  calculation  predicts  Rxx  ,Ryy  and  Rxy  are  significantly 
larger  under  cavitating  conditions  than  under  non-cavitating  conditions  due  to  flow 
reattachment.  It  should  be  mentioned  the  k  —  u  model  does  not  include  any  pro¬ 
duction  terms  due  to  bubble  collapse.  The  only  mechanism  to  capture  turbulence 
generation  for  the  turbulence  model  is  the  occurrence  of  a  strong  strain  rate.  How¬ 
ever,  due  to  the  presence  of  cavitation  the  pseudo  density  model  produces  a  larger 
strain  rate  flow  field.  It  has  been  observed  that  the  collapse  of  bubbles  is  a  source  of 
vortex  generation  (Gopalan  and  Katz,  2000).  Here,  we  can  see  a  good  coordination 
of  the  pseudo  density  model  and  the  turbulence  model  in  terms  of  the  prediction 
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Figure  4.18  (a)  Pseudo  density  contours,  (b)  TKE  (turbulent  kinetic  energy) 
contours,  (c)  eddy  viscosity  contours.  Lj D  =  9.375,  AP  =  28^51 
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Figure  4.19  (a)  Pseudo  density  contours,  (b)  TKE  (turbulent  kinetic  energy) 
contours,  (c)  eddy  viscosity  contours.  LjD  —  9.375,  AP  =  32psi 
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Figure  4.20  (a)  Pseudo  density  contours,  (b)  TKE  (turbulent  kinetic  energy) 
contours,  (c)  eddy  viscosity  contours.  Lj D  =  9.375,  AP  =  36psi 


Figure  4.21  Turbulence  intensities  at  the  exit  under  different  operating  conditions 

{LjD  =  9.375) 


Figure  4.22  Turbulence  intensities  and  mass-averaged  Reynolds  stress  comparison 
with  experimental  results  at  different  axial  locations,  with  and  without  cavitation 

(turbulent  calculation) 
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Figure  4.23  Mean  axial  velocity  comparison  with  experimental  results  at  different 
axial  locations,  with  and  without  cavitation 

of  turbulence  quantities  in  spite  of  the  fact  that  turbulence  production  by  bubble 
collapse  is  not  explicitly  accounted  for. 

Fig.  4.23  shows  the  computed  and  measured  mean  axial  velocity  profiles  at  dif¬ 
ferent  streamwise  locations.  The  velocity  is  normalized  by  the  average  mean 
velocity  at  the  exit.  Both  the  model  and  measurement  have  near  identical  velocity 
profile  for  both  cavitating  flow  and  non-cavitating  flow  at  all  locations  except  inside 
the  cavitation  region(a:  =  1).  However,  the  model  underpredicts  the  potential  core 
velocity  inside  the  slot(a:  =  0, 1,2,3),  especially  for  the  flow  over  the  cavity(a;  =  1). 
The  reason  for  that  might  be  of  the  assumption  of  no  slip  between  phases  employed 
in  the  pseudo  density  model.  In  reality,  the  fluid  over  a  bubble  cavity  tends  to  move 
fa.ster. 

We  also  used  the  laminar  code  to  compute  the  Reynolds  stresses  from  the  unsteady 
solution.  In  order  to  do  that,  the  definition  of  the  Favre  average  of  velocity  was 
utilized. 

ff  ^ 

u  =  u  —  u 


(4.18) 
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Hence,  we  can  express  the  mass-averaged  Reynolds  normal  stress  as 


pu"u"  =  p{u  —  uy 

=  p{v?'  —  2uu  u^) 
=  pu^  —  2upu  -|-  pu^ 
=  pu^  —  2'pv?  d- 
=  pu^  —  'pu^ 


(4.19) 


Similarly,  we  have 

pu'v"  =  puF  —  puv 


n  n  9  - 2 

pv  V  ~  pv^  —  pv 


(4.20) 

(4.21) 


To  complete  the  calculation  of  the  three  components  of  Reynolds  stress  tensor  shown 
above,  the  following  time-averaged  and  mass-averaged  variables  need  to  be  determined 
first  as  follows 
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The  instantaneous  flow  fields  p,  u,  v  from  laminar  calculation  were  stored  to  disk 
during  each  time  step,  and  N  is  the  total  sample  number  of  each  flow  variable  stored. 
Fig.  4.24  shows  the  comparison  of  Reynolds  stresses  from  laminar  calculation  and 
Ruiz  and  He’s  measurements.  As  seen  from  the  figure,  the  disagreement  between  the 
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two  sets  of  data  is  much  bigger,  compared  with  the  results  from  turbulent  calculation. 
Although  the  laminar  calculation  generates  an  unsteady  solution,  it  is  not  capable 
of  capturing  the  very  small  length  scale  turbulence.  On  the  contrary,  regarding  the 
prediction  of  turbulence  intensity  the  turbulence  model  does  a  better  job  on  the  same 
grid  resolution. 

4.4  Flow  through  a  2-D  Nozzle 

Gopalan  and  Katz  (2000)  performed  an  experiment  to  study  the  closure  region 
of  an  attached  cavity  over  a  nozzle  surface.  They  presented  data  on  instantaneous 
and  averaged  velocity,  vorticity  and  turbulence  downstream  of  the  cavitation  region 
when  cavitation  is  very  slight.  The  two  dimensional  code  is  utilized  to  simulate  the 
flow  over  a  nozzle  surface  whose  geometry  is  consistent  with  the  one  in  Gopalan  and 
Katz’s  experiment.  The  cross  section  of  the  water  tunnel  in  the  experiment  is  6.35  x 
5.08  cm^.  Strictly  speaking,  we  should  consider  this  flow  field  as  three-dimensional. 
But  in  order  to  use  our  existing  turbulent  code  to  simulate  this  flow,  we  simplify  this 
problem  to  a  two-dimensional  flow  as  Gopalan  and  Katz  did  in  the  measurement. 

The  computational  mesh  used  is  shown  in  Fig.  4.25.  The  reference  dimension  in 
this  case  is  the  length  of  the  curved  nozzle  surface,  L,  as  shown  in  Fig.  4.25.  A  non- 
uniform  mesh  was  generated,  with  dense  grid  points  near  upper  and  lower  walls  and 
within  the  throat  section  of  the  nozzle.  For  clarity,  only  part  of  whole  computational 
grid  is  shown  in  the  plot,  and  every  other  index  in  both  normal  and  axial  direction  of 
the  grid  is  skipped.  To  ensure  the  accuracy  of  the  calculation  the  inflow  boundary  is 
put  at  a  distance  of  L  upstream  of  the  beginning  point  of  the  nozzle  curve,  and  the 
outflow  boundary  was  put  at  a  distance  of  2L  downstream  of  the  end  point  of  the 
nozzle  curve.  A  grid  refinement  study  similar  to  that  in  section  4.3.2  was  undertaken, 
resulting  in  a  grid  of  140x50  being  adequate  to  resolve  the  flow  field.  The  velocity 
profile  at  xjL  =  0.3  was  chosen  to  test  the  convergence  of  the  solution.  The  result  is 
shown  in  Fig.  4.26. 


Figure  4.24  Turbulence  intensities  and  mass-averaged  Reynolds  stress  comparison 
with  experimental  results  at  different  axial  locations,  with  and  without  cavitation 

(laminar  calculation) 
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Figure  4.25  Computational  mesh  for  the  flow  through  a  nozzle 

At  the  inflow  boundary  uniform  velocity  was  specified.  The  other  boundary  condi¬ 
tions  are  similar  to  that  of  the  flow  through  an  injector  passage.  The  inflow  turbulence 
level  k  is  chosen  as  0.1%,  a  value  close  to  that  in  reality  (Gopalan  and  Katz).  It  is 
found  that  the  solution  of  this  flow  is  dependent  on  the  inflow  value  of  cu.  If  the  inflow 
uj  value  is  set  as  such  that  the  inflow  eddy  viscosity  is  over  thousand  or  more,  the 
unsteadiness  in  cavitation  would  disappear  completely.  To  have  a  better  match  with 
the  measured  cavitation  extent,  the  inflow  value  of  oj  is  chosen  as  such  that  the  bulk 
value  of  the  inflow  eddy  viscosity  is  unity,  a  setup  similar  to  that  of  the  orifice  flow. 
However,  unlike  the  orifice  flow,  under  this  inflow  conditions  an  unsteady  solution  is 
obtained  for  this  flow. 


4.4.1  Pressure  distribution 


First  of  all,  the  pressure  distribution  over  the  nozzle  surface  was  computed  through 
single  phase  runs  and  compared  with  the  computed  results  of  Gopalan  and  Katz 
who  utilized  the  commercial  CFD  code  Fluent,  with  RNG  k  —  e  turbulence  model 
treatment,  to  carry  out  their  computation.  Here,  the  pressure  coefficient  Cp  was 
calculated  in  the  comparison.  It  is  defined  as 


Cp  = 


P  —  Pin 


(4.23) 
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Figure  4.26  Velocity  profiles  at  xj L  —  0.3  for  various  grids 

where  Pin  is  the  inlet  pressure  and  Vin  is  the  velocity  at  the  inlet.  As  specified  by 
the  geometry  of  the  nozzle  surface,  a  particle  of  fluid  will  experience  the  following 
four  typical  phases  (presumably  there  is  no  flow  separation)  as  it  passes  over  the 
surface:  an  increase  in  pressure  near  the  leading  edge,  sharp  decrease  in  pressure 
to  the  minimum  value,  quick  recovery  to  pressure  lower  than  outlet  pressure,  and 
gradual  recovery  to  pressure  close  to  outlet  pressure. 

Fig.  4.27  shows  the  evolution  of  the  Cp  distribution  from  the  k  —  u  model  cal¬ 
culation  and  the  results  of  Gopalan  and  Katz.  As  seen  from  the  plot.  Fluent's  RNG 
calculation  reproduces  all  the  four  phases  mentioned  above.  Although,  in  the  early 
stage  of  the  calculation  (t=0.5)  the  k  —  u  model  generates  a  nearly  identical  pressure 
distribution  as  Fluent  does,  after  a  very  short  time,  the  recovery  of  pressure  down¬ 
stream  oi  xfL  =  0.5,  where  a  strong  adverse  pressure  gradient  is  present,  is  destroyed. 
The  reason  for  that  is  boundary  layer  separation.  We  may  refer  to  figures  4.28  and 
4.29  which  show  velocity  vector  plots  at  t=0.5  and  t=3,  respectively.  They  indicate 
that  there  is  no  any  flow  separation  at  t=0.5,  while  a  big  recirculation  zone  present 
near  the  trailing  edge  at  t=3,  which  destroys  the  pressure  recovery.  It  is  clear  that  the 
k  —  to  model  and  Fluent's  RNG  model  generate  different  pressure  distributions.  We 
would  like  to  reproduce  the  pressure  distribution  as  Fluent  does.  For  this  purpose. 
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Figure  4.27  Pressure  distribution  from  k  —  u  calculation 

we  developed  our  own  RNG  code,  using  the  model  formulation  in  Wilcox(1998),  and 
made  a  run  to  compute  the  pressure  distribution.  The  results  due  to  the  RNG  model 
is  shown  in  Fig.  4.30.  Our  implementation  of  the  RNG  model  has  a  behavior  similar 
to  that  of  the  k  —  u>  model  does,  except  that  it  produces  a  pressure  recovery  closer  to 
Gopalan’s  result.  But  there  still  exists  discrepancy  between  the  two  sets  of  data.  The 
reason  for  that  is  not  clear.  One  possibility  is  that  the  turbulence  model  implemented 
in  Fluent  is  proprietary  and  ,  therefore,  the  results  can  not  be  matched  exactly. 

Another  issue  is  that  Gopalan  and  Katz  did  not  show  a  comparison  of  experimental 
result  and  the  Fluent  result.  Therefore,  it  is  not  certain  that  the  Fluent  results  is 
correct.  Further  more.  Our  computed  results  can  be  partially  validated  by  Gopalan 
and  Katz’s  experiment  in  the  following  aspects.  First,  the  measurement  shows  that 
the  ratio  of  the  maximum  velocity  to  the  inlet  velocity  is  2.4,  whereas  this  ratio  for 
our  computed  velocity  is  2.6.  Second,  there  are  no  boundary  layer  separation  in  the 
high  adverse  pressure  gradient  region  downstream  of  the  minimum  pressure  point  in 
both  the  k  —  u  calculation  and  the  experiment.  However,  given  the  comparison  of 
our  cavitating  flow  calculation  with  experiment,  it  is  clear  that  our  turbulence  model 
predicts  too  large  a  separation  near  the  trailing  edge  of  the  nozzle. 


I 


Figure  4.28  Velocity  vector  with  no  flow  separation,  at  i  =  0.5,  the  k  —  up  model 

calculation 


x/L 

Figure  4.29  Velocity  vector  with  flow  separation  at  i  =  3,  the  k  —  oj  model 

calculation 
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Figure  4.30  Pressure  distribution  from  RNG  model  calculation 


4.4.2  Cavitation  extent 


For  this  flow,  the  cavitation  number  is  defined  as 


K  = 


(4.24) 


Gopalan  and  Katz  observed  that  as  cavitation  number  is  reduced  slightly  below  the 
inception  level,  a  cavity  with  limited  extent  occurs  just  downstream  of  the  minimum 
pressure  point,  with  a  “blunt,  glossy”  leading  edge  and  hairpin-like  structure  down¬ 
stream  of  the  closure  region.  If  the  pressure  level  of  the  system  is  reduced  further 
(decrease  in  cavitation  number),  cavitation  will  develop  into  an  advanced  stage  with 
a  massive  bubble  cloud  shed  downstream. 

In  1995,  to  validate  the  pseudo  density  model  Chen  and  Heister  made  extensive 
comparison  of  cavitation  region  with  experimental  data  for  flows  over  a  conic  head 
and  an  ogival  head  which  are  analogous  to  this  flow,  and  they  got  satisfactory  matches 
with  the  measurements.  To  further  validate  the  turbulent  code,  we  run  it  for  the  flow 
condition  of  KT  =  2.5  under  which  the  flow  reaches  an  advanced  cavitation  stage.  The 
computed  density  contours  for  five  numerical  snapshots  in  time  are  compared  with 
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Gopalan  and  Katz’s  experimental  photographic  shots  of  the  cavitation  region  at  dif¬ 
ferent  instances,  and  the  results  are  shown  in  Fig.  4.31.  In  the  plot,  the  experimental 
snap  shots  have  a  light  region  that  denotes  the  vapor  cavity,  on  the  contrary  the 
computed  density  contours  have  a  dark  region  which  denotes  the  cavity  region  with 
the  outmost  contour  corresponding  to  p  =  0.98.  The  comparison  shows  an  overall 
qualitatively  good  match  between  the  measurement  and  computation.  However,  as 
the  cavity  extends  longer  downstream  (the  two  frames  at  the  bottom  of  the  plot),  the 
experiment  shows  the  bubbly  cloud  region  spreads  wider  both  in  axial  and  normal 
direction  as  compared  with  the  computed  pseudo  density  contours.  The  numerical 
results  seem  to  over  collapse  the  bubble  cloud.  The  reason  for  that  might  be  that  in 
the  wake  of  the  cavitation  region  the  effect  of  vortex  stretching  and  shedding  is  dom¬ 
inant  and  our  calculation  is  two-dimensional  which  does  not  allow  vortex  stretching 
to  occur. 

We  also  performed  a  case  study  of  K  =  4.69  which  is  slightly  lower  than  the 
inception  cavitation  index,  with  inlet  velocity  =  5.2ml s  and  Reynolds  number 
Re  =  1.58  X  10®.  Under  this  situation  an  attached  cavity,  whose  length  varies  slightly 
in  time  but  with  no  massive  detached  cavity  patches,  exists  just  downstream  of  the 
minimum  pressure  point.  To  compare  with  experimental  measurement  for  this  case, 
the  time  average  procedure  similar  to  that  of  section  4.3.6  was  performed  on  the  in¬ 
stantaneous  flow  variables  of  density,  and  velocities.  Thus  all  the  following  analysis 
are  based  on  the  the  averaged  flow  fields.  The  computed  cavity  region  was  com¬ 
pared  with  Gopalan  and  Katz’s  measurement,  and  the  results  are  shown  in  Fig.  4.32. 
The  inset  shows  the  global  position  and  overall  extent  of  the  cavity  which  is  pretty 
small  and  slender.  As  seen  from  the  plot,  the  predicted  cavitation  region  is  in  good 
agreement  with  the  measured  data. 

The  detailed  streamline  near  the  cavity  region  was  shown  in  Fig.  4.33.  The  con¬ 
tours  of  the  pseudo  density  with  p  =  0.6  at  the  outmost  is  shown  in  the  plot.  Because 
of  the  presence  of  the  cavity  the  flow  is  pushed  away  from  the  wall,  and  pulled  back 
at  the  rear  end  of  the  cavity  due  to  the  recovery  of  pressure  there.  A  typical  flow 
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around  the  rear  portion  of 
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structure  around  such  a  cavity  would  be  that  the  flow  separates  from  the  wall  at  the 
leading  edge  of  the  cavity,  passes  over  the  cavity  body,  and  reattaches  to  the  wall  at 
the  trailing  edge  of  the  cavity.  In  other  words,  the  presence  of  cavity  serves  as  an 
extension  of  the  geometry  of  the  surface  which  further  increases  the  adverse  pressure 
gradient  that  might  not  be  strong  enough  to  cause  the  flow  separation  in  a  single 
phase  flow.  Many  researchers  regard  flow  reattachment  as  the  reason  for  the  exis¬ 
tence  of  the  reentrant  flow  under  the  cavity.  Because  the  cavity  is  small  and  slender 
in  this  case,  the  additional  cavity-induced  adverse  pressure  gradient  is  slight.  In  our 
calculation  we  did  not  see  a  reentrant  structure  existing  beneath  the  vapor  region. 
This  is  consistent  with  Gopalan  and  Katz’s  observation  in  their  experiment. 

Fig.  4.34  shows  the  distribution  of  the  Reynolds  normal  stresses  pu"u" ,  pv"v"  and 
shear  stress  pu"v" .  As  seen  from  the  plot,  the  computation  predicts  high  turbulence 
levels  close  to  the  cavity  and  in  the  downstream  region.  As  mentioned  in  the  previous 
section,  this  is  because  of  the  high  strain  rate  in  those  regions  which  is  captured 
by  the  model.  Although,  Gopalan  and  Katz  observed  the  same  trend,  there  exists 
big  discrepancy  in  the  detailed  structure  of  the  contours  between  the  measured  and 
computed  results  in  terms  of  the  manner  and  pattern  of  the  distribution  of  the  stresses. 
The  measured  high  turbulence  level  region  is  closer  to  the  wall  compared  with  the 
computed  results.  Fig.  4.35  shows  the  computed  and  measured  displacement  and 
momentum  thicknesses  of  the  boundary  layer  in  the  closure  region.  Gompared  with 
the  measurement,  the  computation  generate  too  much  displacement  and  momentum 
deficit  in  the  boundary  layer.  It  overpredicts  the  thicknesses  more  than  two  times 
larger  than  the  measured  ones.  Furthermore,  the  typical  form  factor  of  a  turbulent 
boundary  layer  is  about  1.3,  while  the  computed  value  of  this  variable  is  about  1.43. 
The  measured  form  factor  is  about  1.2  at  x/L  =  0.25  —  0.35,  and  it  approaches  1 
further  downstream.  The  cause  of  the  discrepancies  might  be  due  to  the  inherent 
drawback  in  turbulence  modeling,  the  averaging  simplification  which  throws  away  all 
the  instantaneous  information  of  the  flow  field.  One  other  reason  is  because  we  are 


Figure  4.34  Contours  of  Reynolds  stresses  at  the  rear  portion  of  an  attached  cavity: 

(a)/t)u"u",  {h)pu"v" ,  (c)  pv"v" 
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using  a  traditional  single  phase  turbulence  model  which  does  not  take  into  account 
the  effect  of  phase  changes. 


Figure  4.35  Displacement  and  momentum  boundary  layer  thickness  and  comparison 

with  experimental  measurement 


4.5  Conclusions 

The  homogeneous  pseudo  density  model  and  the  k  —  u  turbulence  model  were 
used  to  approximate  cavitation  two-phase  flows.  Simulations  of  the  cavitating  flow 
in  an  injector  slot/oriflce  were  undertaken.  The  result  is  found  sensitive  to  the  inflow 
condition  of  u.  The  inflow  condition  is  set  as  such  that  the  cavitation  length  does 
not  vary  with  the  change  of  this  variable  at  the  inflow  boundary,  and  a  steady  so¬ 
lution  is  obtained.  Turbulence  model  shows  slight  improvement  in  prediction  of  the 
discharge  coefficient.  The  prediction  of  cavitation  length  has  various  performance  for 
different  length/diameter  ratio  slot,  and  worst  agreement  with  experiment  result  is 
obtained  for  the  shortest  slot.  The  computed  mean  velocity  and  Reynolds  stresses 
profiles  show  good  agreement  with  the  measured  results.  Computation  also  indi¬ 
cates  the  trend  of  cavitation  increasing  the  displacement  and  momentum  thicknesses, 
turbulence  intensities  which  is  in  accord  with  experimental  observations. 
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Simulations  of  the  cavitating  flow  through  a  two-dimensional  nozzle  were  also 
performed.  The  computed  pressure  distribution  from  the  k  —  u  model  calculation 
dose  not  agree  with  the  Fluent's  results  due  to  predicting  too  large  a  separation. 
However  no  experimental  results  are  available  to  validate  these  two  calculations.  The 
overall  matches  of  the  computed  cavitation  region  to  the  experimental  photographic 
shot  of  advanced  cavitation  and  measurement  of  attached  cavitation  are  very  well. 
The  computation  also  indicates  high  turbulence  intensity  develops  around  an  attached 
cavity.  However  the  pattern  of  the  distribution  of  the  Reynolds  stresses  does  not  agree 
with  the  measured  ones.  In  addition,  the  models  overpredicts  the  momentum  and 
displacement  deficit  of  the  boundary  layer  in  the  closure  region. 


9  Appendix  C  -  3D  Cavitating  Flow  Simulations 


Xu,  C.,  “Simulation  of  Orifice  Internal  Flows  Including  Cavitation  and  Turbulence”  Ph. 
D.  Thesis,  School  of  Aeronautics  and  Astronautics,  Purdue  University,  West  Lafayette,  IN, 
August  2001,  pp.  82  -  100.  1996. 
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5.  3D  CAVITATING  FLOW  SIMULATIONS 

5.1  Introduction 

Orifice  flows  driven  by  the  cross  flow  in  a  manifold  occur  in  many  atomization 
devices,  such  as  the  injector  of  a  diesel  engine  and  a  liquid  rocket  engine.  The  three- 
dimensional  nature  of  this  particular  flow  and  the  potential  of  it  easily  becoming  cav- 
itated  make  numerical  modeling  of  it  a  great  challenge.  In  2000,  Bunnell  and  Heister 
investigated  the  flow  inside  this  geometry  by  solving  the  fully  three-dimensional,  un¬ 
steady,  two-phase  Naiver-Stokes  equations,  with  cavitation  being  treated  by  the  the 
homogeneous  fluid  model.  They  investigated  the  effect  of  cavitation  on  orifice  mass 
flow,  and  found  cavitation  acts  as  a  slipstream  which  tends  to  decrease  the  discharge 
velocity  at  the  exit,  and  thus  the  rnass  flow  rate.  Experimentally,  Strakey  and  Talley 
(1999)  have  investigated  the  effect  of  manifold  cross  flow  on  the  discharge  coefficient 
of  an  orifice.  They  found  that  the  discharge  coefficient  closely  correlated  with  the 
cross  flow  velocity.  The  discharge  coefficient  can  be  decreased  as  much  as  50%  as  the 
cross  flow  velocity  is  increased  beyond  a  certain  value. 

In  this  chapter,  the  fully  three-dimensional  two-phase  Naiver-Stokes  solver  devel¬ 
oped  by  Bunnell  (1999)  is  utilized  to  simulate  an  orifice  flow  driven  by  a  manifold 
cross  flow.  The  effects  of  the  cross  flow  velocity  on  cavitation  length  and  discharge 
coefficient  are  investigated.  The  characteristics  of  this  internal  flow  field  are  also  an¬ 
alyzed  based  on  the  results  from  the  calculation.  The  computed  discharge  coefficient 
is  compared  with  Strakey  and  Talley’s  measurement.  The  details  of  the  development 
and  numerics  of  the  model  are  described  by  Bunnell  (1999).  Here,  the  results  of  the 
simulation  are  provided  in  the  following  section. 


Figure  5.1  3-D  manifold  cross  flow  model  and  key  parameters,  from  Bunnell  (1999) 
5.2  Case  Descriptions 

Fig.  5.1  shows  the  schematic  representation  of  the  three-dimensional  manifold 
cross  flow.  In  the  figure,  Vi  is  the  cross  flow  velocity,  and  V2  is  ideal  discharge  velocity 
at  the  exit,  calculated  from  Bernoulli’s  equation.  In  this  study,  eight  simulations  were 
conducted  on  an  orifice  with  diameter  D  =  2.03mm  and  aspect  ratio  LjD  =  b.  In 
all  the  cases,  the  upstream  pressure  is  held  as  Pi  =  O.Q9MPa,  a  value  which  was 
used  in  Strakey  and  Talley’s  experiment.  By  varying  the  back  pressure  of  the  orifice, 
hence  the  cavitation  number,  and  the  cross  flow  velocity,  we  get  the  various  operating 
conditions  used  in  the  simulations.  They  are  shown  in  Table  5.1  (Note  cases  4,  5, 
and  6  were  simulated  by  Bunnell.  For  completeness,  we  list  all  of  them  in  the  table). 

The  operating  conditions  are  consistent  with  those  used  by  Strakey  and  Talley  in 
their  experiment.  Therefore  it  is  possible  for  us  to  make  comparison  of  the  predicted 
discharge  coefficient  with  the  measurements.  All  the  simulations  were  performed 
using  a  mesh  of  120x42x42  as  recommended  by  Bunnell.  In  order  to  save  computing 
time,  for  each  simulation  the  first  30,000  time  steps  were  carried  out  as  single  phase 
flow  calculation,  and  after  that  it  was  run  as  two  phase  flow  calculation.  A  typical 
run  takes  3-4  weeks  on  two  Pentium  650HZ  PCs.  Because  it  is  very  time-consuming, 
only  5  simulations  were  performed  in  this  study.  Those  data  and  the  three  cases  run 
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Table  5.1  Operating  conditions  simulated 


Case 

I/i(m/s) 

K 

LcfL 

1 

6.0 

1.2 

0.56 

2 

6.0 

1.6 

0.04 

3 

6.0 

6.0 

0 

4* 

8.9 

1.2 

0.72 

5* 

8.9 

1.8 

0.02 

6* 

8.9 

6.0 

0 

7 

12.1 

1.2 

1 

8 

15.1 

1.2 

1 

*  from  Bunnell  (1999) 


by  Bunnell,  complete  the  comparison  with  Strakey  and  Talley’s  measurement  for  one 
orifice  geometry  (L/D  =  5). 

5.3  Results 

The  cavitation  length  history  for  cases  1  and  2  are  shown  in  Fig.  5.2.  Here,  Lc 
is  the  length  of  the  cavitation  zone  as  defined  as  the  most  downstream  point  where 
pseudo  density  is  0.98.  As  found  before,  cavitation  number  is  an  important  parameter 
in  determining  whether  or  not  the  flow  becomes  cavitated.  If  cavitation  number  is 
low,  ie.  the  discharge  pressure  level  is  close  to  the  vapor  pressure,  the  cavitation 
inside  the  orifice  becomes  prominent.  The  computed  averaged  cavitation  length  for 
each  cases  is  shown  in  Table  5.1.  In  this  series  of  simulations,  when  K  =  1.2,  for  all 
cases  where  cross  flow  velocity  ranges  from  6.0mls  to  15.1m/s,  a  highly  developed 
cavity  occurs  inside  the  orifice.  As  cavitation  number  increases  to  1.6  (Vi  =  6.0m/s) 
or  1.8  {Vi  —  S.9m/s),  cavitation  reduces  dramatically  in  both  cases.  With  a  further 
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increase  in  this  number  to  6.0,  the  flow  remains  non-cavitating  completely  for  any  of 
the  values  of  cross  flow  velocity  simulated  in  this  study. 


t 

Figure  5.2  Cavitation  length  history,  Vi  =  6.0m/ s 

Fig.  5.3  shows  the  cavitation  length  history  for  various  cross  flow  velocities  when 
K  —  1.2.  As  seen  from  this  plot,  the  other  point  is  evident  that  cavitation  becomes 
more  significant  as  the  flow  moves  faster  across  the  orifice  in  the  manifold.  This 
happens  because  with  increased  cross  flow  velocity  the  separated  region  at  the  leeward 
corner,  where  cavity  develops,  becomes  larger.  This  trend  agrees  with  the  observations 
of  Strakey  and  Talley.  However,  the  dominant  factor  that  determines  the  extent 
of  cavity  is  still  the  cavitation  number.  This  point  can  be  verified  by  comparing 
cavitation  length  of  cases  2  and  5.  Case  2  corresponds  to  a  larger  cavity  because 
of  the  lower  cavitation  number,  even  though  its  cross  flow  speed  is  less  than  that 
of  case  5.  Strakey  and  Tally  found  that  the  cavitation  inception  index  is  about  1.8. 
Due  to  the  high  expense  of  this  computation,  we  can  not  afford  to  determine  the 
the  exact  inception  cavitation  number,  but  the  computed  cavity  is  very  small  when 
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Figure  5.3  Cavitation  length  history  for  various  cross  flow  velocity,  when  K  =  1.2 

cavitation  number  used  in  the  simulations  (cases  2  and  cases  5)  are  equal  or  close  to 
the  measured  inception  index. 

Fig  5.4  shows  the  history  of  the  discharge  coefficient  for  the  cases  of  Vi  =  Q.OmIs. 
Bunnell  and  Heister  had  a  thorough  analysis  of  the  unsteady  characteristics  of  the 
discharge  coefficient  before.  Here  we  observe  the  same  periodic  behavior  of  the  dis¬ 
charge  coefficient  as  they  did.  As  seen  from  the  plot,  the  oscillation  in  the  discharge 
coefficient  still  exists  for  the  non-cavitating  case  due  to  the  instability  induced  by  the 
vena-contracta.  For  the  cavitating  cases,  after  carefully  examining  the  curves  of  the 
discharge  coefficient  and  that  of  their  corresponding  cavitation  length,  we  can  see  the 
instant  at  which  the  discharge  coefficient  reaches  its  peak  lies  in  the  increasing  phase 
of  the  cavitation  length,  somewhat  slightly  later  than  the  instant  at  which  cavitation 
length  reaches  its  peak.  This  indicates  that  the  increase  in  the  volume  of  the  cav¬ 
ity  pushes  a  larger  amount  of  liquid  out  of  the  exit,  on  the  other  hand  during  the 
shrinking  phase  of  the  cavity  more  liquid  is  refrained  inside  the  orifice. 
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Figure  5.5  shows  the  comparison  of  the  computed  average  discharge  coefficient 
and  the  experimentally  measured  discharge  coefficient.  Note  the  average  value  is  ob¬ 
tained  by  sampling  the  data  after  it  reaches  a  steady  oscillation  state  (For  the  case  of 
K”  =  1.2,  we  should  run  the  code  longer.  But  as  seen  from  the  plot,  a  nearly  steady 
oscillation  state  is  obtained  for  this- case.  The  average  discharge  coefficient  for  this 
case  is  obtained  by  sampling  the  data  between  t— 116-150).  Although  the  number 
of  computed  data  is  fewer  than  that  of  the  experimental  results  because  of  the  high 
expense  of  the  computation,  two  patterns  can  still  be  observed  from  the  plot.  First, 
for  a  given  cross  flow  velocity,  as  cavitation  number  increases  the  discharge  coeffi¬ 
cient  increases  at  first.  After  reaching  the  maximum  value,  discharge  coefficient  then 
begins  to  drop  off.  The  maximum  value  in  discharge  coefficient  occurs  when  the  cav¬ 
itation  number  is  about  1.8,  a  value  close  to  the  experimentally  measured  inception 
index.  The  reason  for  that  is  the  following.  The  mass  flow  rate  is  proportional  to 
the  pressure  drop  through  the  orifice  (increasing  with  decreasing  cavitation  number) 
before  the  inception  of  cavitation.  When  the  pressure  loss  is  large  enough  to  cause 
cavitation  to  occur,  the  mass  flow  rate  begins  to  drop  due  to  the  slipstream  effect  of 
the  cavity  as  pointed  out  by  Bunnell  and  Heister.  The  more  significant  cavitation 
becomes,  the  lower  the  mass  flow  rate  is.  The  limited  computed  data  samples  repro¬ 
duce  this  feature  as  the  experimental  counterparts  do,  except  that  they  overpredict 
the  discharge  coefficient  in  general,  especially  when  Vi  =  8.9m/s  and  K  =  6.0.  The 
overprediction  might  be  due  to  the  lack  of  turbulence  in  the  modeling  or  because  the 
grid  is  not  fine  enough  (Bunnell,  1999). 

The  second  trend  observed  in  Fig.  5.5  is  with  regard  to  the  effect  of  cross  flow 
velocity.  Both  the  computation  and  experiment  show  that  as  cross  velocity  increases, 
discharge  coefficient  decreases.  This  phenomenon  can  be  attributed  to  the  increase  in 
velocity  ratio,  the  ratio  of  the  cross  flow  velocity  to  the  ideal  orifice  discharge  velocity 
Vi/V2,  (Strakey  and  Talley,  1999).  With  increasing  velocity  ratio,  less  amount  of  fluid 
passes  through  the  orifice. 
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Figure  5.5  The  effect  of  cross  flow  velocity  on  discharge  coefficient  and  the 
comparison  with  experiment  measurement 


81 


Table  5.2  Computed  average  discharge  coefficient,  its  RMS  oscillation  and 
oscillating  frequency,  with  cross  flow  velocity  Vi  =  6.0m/s 


K 

Cd 

<1 

f  (KHZ) 

1.2 

0.697 

0.00616 

1.50 

1.6 

0.768 

0.00165 

3.69 

6.0 

0.722 

0.00163 

2.00 

Table  5.2  shows  the  computed  average  discharge  coefficient  as  well  as  its  RMS 
oscillation  and  oscillating  frequency  for  the  cases  where  Vi  =  6.0m/s.  The  data  show 
a  pattern  similar  to  what  we  found  for  the  axisymmetric  flows  (Chapter  2),  when 
cavitation  occurs  longer  cavities  tend  to  oscillate  at  lower  frequency.  Also  note  that 
the  variation  of  mass  flow  is  increased  dramatically  for  highly  cavitated  flow  because 
the  oscillation  of  the  size  of  the  bubble  cavity  is  the  primary  cause  of  the  variation. 

To  help  us  understand  how  the  three-dimensional  bubble  cavity  evolves,  the  iso¬ 
surface  on  which  the  pseudo  density  is  0.98  and  the  slice-contour  of  the  pseudo  density 
at  different  instants  in  time  are  shown  in  figures  5.6  through  5.9  and  figures  5.10 
through  5.13,  respectively.  In  figures  5.6  through  5.9,  the  cross  flow  is  from  left  to 
right  and  green  lines  are  streamlines.  In  figures  5.10  through  5.13,  the  location  and 
orientation  of  the  contour  planes  are  indicated  in  the  three-dimensional  inset  in  the 
figures,  and  for  the  angular  plane  the  flow  is  from  right  to  left. 

The  four  instants  indicated  in  the  plots  lie  within  one  cycle  of  the  quasi-periodic 
oscillation  of  the  cavity  (refer  to  Fig.  5.2).  As  seen  from  figures  5.10  through  5.13,  the 
shape  of  the  forepart  of  the  cavity  is  nearly  unchanged  in  time  in  both  circumferential 
and  radial  directions.  However,  in  the  wake  region,  the  shape  and  extent  of  the  rear 
part  of  the  cavity  varies  considerably  with  time.  In  the  growth  phase  of  the  cavitation 
length  (t  =  122,  124),  cylindrical  like  structures  exist  at  the  rear  part  of  the  cavity. 
This  is  similar  to  what  was  found  by  researchers  (Kubota,  et  ah;  1992),  based  on 
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the  observation  of  the  cavitating  flow  over  a  hydrofoil,  that  the  breakup  of  sheet 
cavitation  generates  large  cylindrically  shaped  bubble  clouds  shed  downstream. 

Also  note  the  three  stream  lines,  which  begins  at  flxed  locations.  The  changes 
of  these  three  specific  streamlines  highlight  the  variation  of  the  flow  structures  of 
this  cavitating  flow  during  a  cycle  of  the  growth  and  shrinking  of  the  cavity.  The 
streamline  on  the  windward  side  is  straight  and  nearly  unchanged  with  time,  however 
the  one  on  the  leeward  side  is  highly  sinuous.  The  sinuousness  or  twisting  of  the 
streamline  on  the  leeward  side  indicates  large  eddies  are  generated  downstream  of 
bubble  cavity. 

Soteriou  et  al.  (2001)  have  found  strong  vortex  structures  exist  inside  the  orifice 
when  cavitation  occurs  based  on  their  observation  on  a  large  scale  nozzle.  Vortex 
interaction  can  also  be  observed  under  non-cavitating  conditions.  The  vortices  inside 
the  nozzle  sometimes  intertwine  with  each  other.  Fig.  5.14  shows  the  photographic 
shot  of  a  vortex  intertwining  structure  inside  an  orifice,  and  Fig.  5.15  shows  the  vortex 
structure  of  a  diesel  spray.  Contrary  to  the  traditional  theory  that  the  aerodynamic 
interfacial  shear  is  the  dominant  force  that  causes  the  breakup  of  the  fuel  injected  from 
the  orifice,  Soteriou  et  al.  indicated  that  the  primary  factors  inducing  atomization 
are  associated  with  the  characteristics  of  the  internal  flow  inside  the  nozzle. 


Figure  5.14  An  experimental  photographic  snapshot  of  vortex  structure  inside  an 
orifice.  [From  Soteriou  et  al.  (2001),  used  by  permission] 
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Figure  5.15  An  experimental  photographic  snapshot  of  vortex  structure  of  diesel 
sprays,  on  the  right  of  the  plot  the  dominant  structure  is  outlined.  [From  Soteriou 

et  al.  (2001),  used  by  permission] 

Fig.  5.16  and  5.17  show  the  various  computed  vortex  lines  inside  the  orifice  at  one 
instant  in  time  for  the  conditions  with  and  without  cavitation,  respectively.  At  other 
instants,  both  flow  conditions  show  a  vortex  structure  similar  to  what  are  shown  in 
these  two  figures.  The  three  component  of  vorticity  are  in  line  with  streamwise,  cross 
and  vertical  directions.  The  source  that  generates  vorticity  is  the  shear  stress  in  the 
viscous  boundary  layer.  The  vortex  structures  exist  inside  the  orifices  indicates  that 
the  vorticity  generated  from  the  boundary  layer  upstream  of  the  separation  as  well  as 
that  generated  from  the  local  boundary  layer  is  convected  into  the  fluid  in  the  core 
region. 

Note  that  each  vortex  line  in  Fig.  5.16  has  a  corresponding  counterpart  in  Fig.  5.17 
that  starts  at  or  goes  through  the  same  location.  However,  as  seen  from  the  plots, 
the  difference  between  the  vortex  structure  under  cavitating  condition  and  that  under 
non-cavitating  condition  is  significant.  For  non-cavitating  case,  in  the  front  section 
of  the  orifice,  the  cross  or  vertical  components  of  the  vorticity  are  dominant.  The  few 
vortex  lines  shown  in  the  figure  are  nearly  parallel  to  each  other  and  perpendicular 
to  streamwise  direction.  Based  on  these  representing  vortex  lines,  we  believe  vortex 
sheet-like  structures  exist  near  the  wall  region  in  the  front  section  of  the  orifice.  In  the 
rear  section  of  orifice,  the  streamwise  component  of  the  vorticity  start  to  be  dominant. 
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Figure  5.16  Vortex  structure  inside  the  orifice  under  non-cavitating  condition, 
Vi  =  6.0m/s,  K  =  1.2.  Green  lines  are  vortex  lines  outlined  parallel  to  local 

vorticity  vector 


Figure  5.17  Vortex  structure  inside  the  orifice  under  non-cavitating  condition, 

Vi  =  6.0m/5,  K  =  6.  Green  lines  are  vortex  lines  outlined  parallel  to  local  vorticity 

vector 
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Cavitating  K  =  1 .2 


Non-cavitating  K  =  6.0 


Figure  5.18  Exit  plane  velocity  streamlines,  Vi  =  6.0m/s.  The  right  of  circles  is  the 

leeward  side. 

The  vortex  sheet  starts  to  roll  up  and  several  vertex  lines  intertwine  and  form  a  strong 
vortex  structure  on  the  windward  side  of  orifice. 

For  the  flow  with  cavitation,  vortex  interaction  begins  earlier  on  the  leeward  side 
and  two  primary  vortex  structures  are  formed  at  the  exit  of  orifice.  These  vortex 
structures  indicates  that  swirls  with  rotating  axis  in  line  with  streamwise  direction 
exist  at  the  exit.  This  phenomenon  has  been  observed  by  Bunnell  and  Heister,  and 
can  be  evidenced  by  figure  5.18  which  shows  the  streamlines  of  the  cross  and  vertical 
components  of  the  velocity  vector  on  the  exit  plane.  It  is  interesting  to  note  that 
more  than  one  swirls  exist  on  the  exit  plane  for  non-cavitating  condition.  Fig.  5.19 
shows  the  contour  of  the  streamwise  component  of  the  vorticity  vector,  u>x,  at  the  exit 
plane.  Combining  figures  5.18  and  5.19,  we  can  see  as  counter-rotating  swirls  exist 
in  the  upper  part  and  lower  part  of  the  exit  plane  the  sign  of  changes  accordingly. 

Fig.  5.20  shows  the  exit-plane-contour  of  the  magnitude  of  the  vorticity  vector, 
which  is  defined  as 

I  ^  1=  +^y  (^-1) 

where  and  are  the  three  components  of  the  vorticity  vector.  Note  the 

maximum  or  minimum  contour  value  of  the  three  variables  plotted  in  these  two  plots 


Cavitating  K  =  1 .2 


Non-cavitating  K  =  6.0 


Figure  5.19  Contour  plot  of  the  streamwise  component  of  vorticity  vector  on  exit 
plane,  Vi  =  Q.Om/s.  The  right  of  circles  is  the  leeward  side. 


91 


Cavitating  K  s  1 .2 


Non-cavitating  K  =  6.0 
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Figure  5.20  Contour  plot  of  the  magnitude  of  vorticity  vector  on  exit  plane, 

Vi  =  6.0m/5.  The  right  of  circles  is  the  leeward  side. 

occurs  near  the  wall  region.  To  highlight  the  structure  of  the  interior  domain,  the 
contour  minimum  and  maximum  levels  have  been  reduced  such  that  the  interior 
structure  is  obvious.  As  seen  from  Fig.  5.20,  for  cavitating  flow,  the  vorticity  on 
the  leeward  side  is  much  higher  than  that  on  the  windward  side.  On  the  contrary 
this  distribution  seems  even  for  the  non-cavitating  flow,  but  with  considerable  smaller 
magnitude  of  vorticity  vector  on  the  leeward  side.  Also,  flgures  5.18  and  5.20  indicate 
that  at  the  center  of  a  swirl  the  magnitude  of  the  vorticity  vector  does  not  reaches 
its  maximum  value. 

The  contours  for  the  magnitude  of  the  vorticity  vector  on  the  angular  plane  in¬ 
clined  45°  with  respect  to  the  cross  flow  vector  at  an  instant  in  time  under  cavitating 
and  non-cavitating  conditions  are  shown  in  Fig.  5.21.  The  flow  is  from  left  to  right 
in  the  figure.  As  seen  from  this  plot,  the  flow  turns  into  the  orifice  with  very  low- 
level  vorticity  in  the  core  region.  But  after  the  separation,  the  magnitude  of  the 
vorticity  in  the  interior  part  of  the  domain  increased  greatly  for  both  cavitating  and 
non-cavitating  flows.  The  difference  made  by  cavitation  is  that  the  vorticity  at  the 
exit  is  stronger  for  a  cavitating  flow. 
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Figure  5.21  Contour  plot  of  the  magnitude  of  vorticity  vector  on  angular  plane, 

Vi  =  6.0m/s. 

5.4  Conclusions 

A  series  of  numerical  simulations  have  been  performed  to  study  the  flow  inside 
an  orifice  driven  by  a  manifold  cross  flow.  The  effect  of  cavitation  number  and 
cross  flow  velocity  on  cavitation  length  and  mass  flow  rate  has  been  examined.  The 
cavitation  number  is  the  primary  factor  that  indicates  to  what  degree  the  flow  would 
become  cavitated.  However,  the  cross  flow  velocity  also  has  a  considerable  effect  on 
the  extent  of  the  cavitation  region.  For  a  given  cavitation  index  at  which  cavitation 
would  occur,  increasing  cross  flow  velocity  increases  cavitation  and  can  promote  a 
hydraulic  flip  condition.  The  presence  of  cavitation  acts  as  a  slipstream  that  tends 
to  reduce  the  mass  flow  rate  through  the  nozzle,  and  a  more  advanced  cavitation 
situation  corresponds  to  less  mass  flow  rate.  Therefore,  an  increase  in  cross  flow 
velocity  causes  a  decrease  of  the  mass  flow  rate  through  the  orifice. 

The  flow  field  inside  the  orifice  under  both  cavitating  and  non-cavitating  condi¬ 
tions  has  also  been  assessed.  The  three-dimensional  cavity  grows  and  shrinks  quasi- 
periodically  with  cylindrical  structures  formed  at  its  rear  part.  These  structures  can 
be  interpreted  as  bubble  clouds  shed  downstream.  Streamwise  vortex  intertwining 
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structure  has  been  found  for  both  cavitating  and  non-cavitating  flows.  However,  the 
vortex  interaction  tends  to  start  earlier  as  the  flow  begins  its  journey  inside  the  orifice 
with  the  presence  of  cavitation.  The  role  played  by  the  internal  vortex  structure  on 
spray  atomization  is  not  well  understood  since  it  is  both  a  numerical  and  an  experi¬ 
mental  new  discovery.  Recently  some  authors  (Hiroyasu,  2000;  Soteriou  et  al.,  2001; 
Tamaki  et  al.,  1998)  have  suggested  that  cavitation  inside  the  orifice  is  an  important 
factor  that  promotes  atomization.  This  might  be  attributed  to  the  difference  made 
by  cavitation  on  the  internal  flow  structure  and  the  magnitude  of  the  vorticity  vector 
at  the  exit  plane.  The  simulation  shows  that  there  is  a  stronger  swirl  on  the  leeward 
side  for  cavitating  flow.  Once  the  fuel  flows  out  of  the  orifice,  without  the  constraint 
of  the  wall  boundary  the  stronger  velocity  components  in  cross  and  vertical  directions 
tend  to  cause  the  jet  to  break  up  earlier. 


